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Abstract 


In  the  classical  form  of  the  finite  element  method,  called  the  h 
version,  piecewise  polynomials  of  fixed  degree  p  are  used  and  the  mesh  size  h 
is  decreased  for  accuracy.  In  this  paper,  we  discuss  the  fundamental 
theoretical  ideas  behind  the  relatively  recent  p  version  and  h-p  version.  In 
the  p  version,  a  fixed  mesh  is  used  and  p  is  allowed  to  increase.  The  h-p 
version  combines  both  approaches.  We  describe  and  explain  the  basic 
properties  and  characteristics  of  these  newer  versions,  especially  in  areas 
where  their  behavior  is  significantly  different  from  that  of  the  h  version. 

We  include  simplified  proofs  of  key  concepts  and  provide  computational 
illustrations  of  several  results. 


Subject  Classifications:  AMS  (MOS):  65N30 


ii 


1, 


INTRODUCTION 


A  major  tool  in  computational  mechanics  today  is  the  finite  element 
method  (FEM).  Although  this  method  has  a  long  history  (see  e.g.,  [Od][Wl]) 
the  beginning  of  its  computational  success  can  be  related  to  [TC] ,  [Cl]. 
Originally,  the  main  application  of  the  FEM  was  in  the  area  of  structural 
mechanics  but  today,  its  use  has  expanded  to  various  other  fields  such  as 
fluid  mechanics,  thermal  analysis,  electrical  engineering,  etc. ,  with 
applications  to  both  linear  and  nonlinear,  stationary  and  transient  analysis. 
Consequently,  there  is  a  vast  body  of  published  research  in  the  field  of 
finite  element  theory  and  applications.  In  Figure  1.1,  we  have  presented  data 
showing  the  number  of  papers  related  to  the  FEM.  (See  also  [Ma].) 
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Figure  1.1.  Number  of  FEM  papers  in  MAKEBASE  database  (1990,  1991  are 
not  completely  correct-delayed  literature  acquisition.  Computational  fluid 
mechanics  is  not  included. )  (Courtesy:  Jaroslav  Mackerle,  Linkoping  Inst,  of 
Technology,  Sweden). 
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Most  of  the  papers  represented  in  Figure  1.1  are  related  to  the 
classical  form  of  the  FEM,  called  the  h  version,  in  which  polynomials  of  fixed 
degree  p  are  used  and  the  mesh  is  refined  to  increase  accuracy.  There  are 
several  commercial  finite  element  programs  based  on  the  h  version  of  the  FEM 
which  are  currently  available  (e.g.,  MSC/NASTRAN,  Cosmos/M,  Abaqus,  Aska, 
Adina,  Ansys,  etc.). 

In  contrast,  only  a  few  (less  than  100)  of  the  above  papers  are  related 
to  the  relatively  recent  p  and  h-p  versions.  The  p  version  of  the  FEM  uses  a 
fixed  mesh  but  increases  the  polynomial  degree  p  to  increase  accuracy.  In  the 
h-p  version,  the  two  approaches  of  mesh  refinement  and  degree  enhancement  are 
combined.  (The  h  and  p  versions  can  be  considered  as  special  cases  of  the  h-p 

version.)  The  first  theoretical  papers  on  these  new  versions  ( [BSK] ,  [BD]  ) 
were  published  in  1981.  Recently,  a  few  commercial  and  large  research 
programs  have  become  available  on  these  versions,  such  as  MSC/PROBE,  Applied 
Structure,  PHLEX  and  STRIPE.  In  contrast  to  the  many  books  devoted  to  the  h 
version  FEM  in  engineering  and  mathematics,  the  only  book  addressing  the  p  and 
h-p  versions  is  [SzB]. 

In  this  paper,  we  discuss  the  basic  theoretical  ideas  of  the  p  and  h-p 
versions  of  the  FEM.  Rather  than  an  exhaustive  summary  of  published  work,  our 
goal  here  is  to  describe  only  the  most  salient  (and  practically  useful) 
features  which  are  characteristic  of  the  p  and  h-p  versions.  We  give 
numerous  elaborations  and  examples  of  key  theoretical  results  to  help  explain 
the  underlying  reasons  for  the  behavior  of  these  versions,  and  to  bring  out 
the  differences  that  exist  in  comparison  to  the  classical  h  version.  We  also 
include  Informal  proofs  of  several  Important  theorems  to  make  them  more 
accessible  and  to  help  introduce  the  reader  to  p  version  theoretical 
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techniques  of  approximation  (which  can  be  quite  different  from  standard  h 
version  techniques).  Throughout  this  survey,  we  concentrate  on  problems 
related  to  structural  mechanics,  as  opposed  to  other  applications. 

The  plan  of  our  paper  is  as  follows.  In  Section  2,  we  consider  the 
one-dimensional  case,  which  allows  us  to  present  various  ideas  and  results 
(which  hold  for  higher  dimensions  as  well)  in  a  simplified  setting.  Perhaps 
the  most  important  characteristic  of  the  p  and  h-p  versions,  the  enhanced 
asymptotic  rates  of  convergence  that  are  possible,  is  discussed  at  length  in 
Sections  2.3  -  2.4.  We  show,  by  elementary  arguments,  exactly  what  leads  to 
the  well-documented  phenomenon  of  the  "doubling"  in  the  rate  of  convergence  of 
the  p  version  (over  the  h  version)  in  the  presence  of  singularities.  We  also 
demonstrate  (by  meai,s  of  some  simple  proofs)  the  following  basic  difference: 
the  rate  for  the  h  version  can  never  be  better  than  algebraic  (no  matter  which 
mesh  is  used),  while  the  h-p  version  can  yield  exponential  rates.  In  Sections 
2.5,  2.6,  we  discuss,  respectively,  the  error  in  the  norm  and  the  problem 
of  "pollution",  two  areas  in  which  the  p  version  behaves  differently  from  the 
h  version. 

The  two-  and  three-dimensional  analogs  of  results  in  Section  2  are 
discussed  in  Section  4.  The  refined  approximation  theorems  for  the  p  and  h-p 
versions  presented  in  Section  4  are  extremely  dependent  on  precise  regularity 
results  for  the  solution  (which  are  more  detailed  than  those  required  for  the 
h  version).  Therefore,  Section  3  is  devoted  to  a  survey  of  various  types  of 
regularity  theorems  for  elliptic  problems  that  are  used  in  the  development  of 
such  convergence  results.  In  particular,  regularity  results  in  terms  of 
countably  normed  spaces,  which  are  crucial  in  deriving  exponential  rates  of 
convergence  for  the  h-p  version,  are  presented  in  Section  3.3. 
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An  inherent  property  of  the  p  version  is  its  robustness  in  most 
situations,  a  characteristic  that  can  be  put  to  good  use  in  applications  to 
problems  that  are  susceptible  to  numerical  locking.  We  explain  this 
phenomenon  in  Section  5  and  survey  some  recent  related  results,  using  the 
equations  of  linear  elasticity  as  a  model. 

The  convergence  results  in  Section  2,  4  are  mostly  in  the  energy  norm. 
Often,  instead  of  the  energy,  some  other  functional  of  the  solution  is  of 
interest.  In  Section  6,  we  discuss  extraction  procedures  by  which  these 
functionals  can  be  recovered  from  the  solution.  As  a  result,  the  high  rates 
of  convergence  for  the  h-p  version  in  the  energy  are  preserved  for  these 
quantities  of  interest  as  well. 

In  Section  7,  we  consider  the  use  of  the  p  version  in  plate  modeling. 

We  show  how  a  hierarchy  of  plate  models  may  be  set  up,  which  approximates  the 
three-dimensional  plate  with  increasing  accuracy.  The  implementation  of  such 
a  hierarchy  of  plate  models  can  then  be  effected  in  a  natural  way  by  using  the 
p  version.  This  may  be  used  to  ensure  convergence  to  the  solution  of  the 
three-dimensional  plate,  as  opposed  to  that  of  a  two  dimensional  plate  model 
which  is  what  happens  when  a  fixed  model  (e.g.,  the  Reissner-Mindlin  plate)  is 
discretized. 

Implementational  aspects  of  the  p  and  h-p  version  are  briefly  discussed 
in  Section  8. 

We  have  not  tried  to  compose  an  exhaustive  bibliography  of  publications 
on  the  p  and  h-p  version  (as  could,  for  instance,  be  obtained  using 
MAKEBASE).  Rather,  we  have  provided  only  selected  representative  references, 
for  the  interested  reader  to  obtain  more  information  on  the  topics  discussed. 
In  this  connection,  some  representative  references  on  the  p  and  h-p  versions 
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for  areas  not  discussed  in  this  paper  are  as  follows:  boundary  element 
methods  [PS],  [SS] ,  mixed  methods  [JV],  [SI],  parabolic  problems  [BJ],  free 
boundary  problems  [ZB],  shape  optimization  [Sh] ,  non-linear  problems  [Ng] , 
solution  of  large  systems  [Mn] .  In  addition,  see  [Sz]  for  a  survey  on 
computational  questions  related  to  the  p  and  h-p  version. 

2.  The  h,  p  AND  h-p  VERSIONS  IN  ONE  DIMENSION 

The  one-dimensional  setting  allows  a  simple  and  convenient  vehicle  for 
the  exposition  of  our  basic  ideas  and  results.  The  fundamental  properties  of 
and  differences  between  the  h,  p  and  h-p  versions  are  most  clearly  explained 
and  understood  through  the  one-dimensional  results  and  numerical  illustrations 
in  this  section.  These  also  indicate  what  can  be  expected  in  two  and  three 
dimensions  (discussed  in  succeeding  sections),  for  which  strongly  analogous 
results  hold.  A  detailed  analysis  of  the  1-d  case  has  been  made  in  [GiB], 
from  which  most  of  the  results  in  this  section  have  been  cited. 

2.1.  The  model  problem  and  its  discretization 

Let  us  consider  the  model  problem 

(2.1)  -u"  =  f  on  Q  =  (0,1)  =  I, 

u(0)  =  u(l)  =  0. 

k-2  1 

We  will  consider  the  general  case  when  it  is  only  known  that  f  e  H  (Jl)  , 

k  £  1,  as  well  as  the  particular  case  when  f  is  such  that  the  exact 
solution  of  (2.1)  is 

^Here  and  elsewhere,  W^’^(S2),  H^(Q),  9^(0)  will  denote  the  usual 
Sobolev  spaces. 
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(2.2)  u(x)  =  x“  -X,  a  >  2  . 

It  will  be  seen  shortly  that  much  more  refined  convergence  results  exist 

for  the  particular  case  (2.2)  than  the  general  case.  The  form  (2.2)  is 

obviously  singular  at  the  endpoint  x  =  0.  The  reason  we  choose  it  is  that  it 

serves  as  a  model  for  two  (and  higher)  dimensions  where  analogous 

singularities  exist  at  the  corners  of  domains  which  are  curvilinear  polygons 

(see  Section  3).  The  restriction  “  ^  ^  is  made  to  ensure  that  u  e  H^(£l). 

In  general,  the  regularity  of  the  solution  u  is  determined  by  the 

k-2  k 

regularity  of  f.  Obviously,  if  f  e  H  (Q),  then  u  e  H  (0).  An  analogous 
but  more  complicated  relation  between  the  regularity  of  the  input  data  and  of 
the  solution  occurs  in  two  and  three  dimensions  and  will  be  discussed  in 
Section  3. 

As  a  prelude  to  using  the  finite  element  method,  the  problem  (2.1)  is 
cast  into  the  following  variational  form.  Find  u  e  ((^(Q)  such  that  for  any 

V  €  Hj(n), 

(2.3)  B(u,v)  =  r  u' (x)v' (x)dx  =  f(v) 
with 

f(v)  =  r  f(x)v(x)dx. 

The  above  may  now  be  discretized  by  constructing  a  sequence  of  finite¬ 
dimensional  (finite  element)  subspaces  c  H!,(n)  of  dimension  N  = 

0  n 

N(V"),  n  =  1,2,...,  called  the  number  of  degrees  of  freedom.  These  spaces 

consist  of  piecewise  polynomials  on  a  sequence  of  meshes  {7^>,  n  =  1,2 . 

on  Q  where 

=  (0  =  x|?<x^<***<x^  =1). 

0  1  m 

n 
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I.  .  i  rnn,.  nn.,  ,  , 

We  denote  I.  =  [x.,,x.J,  h  .=x.-x.,,j  =  l,...,m,  h  =  max  h 

J  J-1  J  n.j  J  J-1  n  n  n.j 

and  call  Ij  elements.^ 

Let  us  now  consider  three  basic  types  of  meshes,  which  are  analogous  to 
those  used  in  the  two  and  three  dimensional  settings. 

1.  The  guasiuniform  mesh.  Here  is  such  that  we  have  for  all  n 

h 

n 

s  T, 


h  . 
n.j 


where  t  >  0  will  be  called  the  quasiuniformity  constant.  A  special  case  is 
the  uniform  mesh,  where  x  =  1  and 


(2.5) 


n 

X. 

1 


=  —  ,  i  =  0 . m 


n 


2.  The  radical  mesh  with  power  0  >  0  (with  respect  to  x  =  0).  This 
is  defined  by  choosing 


(2.6) 


i  =  1, 


,  m 


(the  uniform  mesh  is  obviously  a  radical  mesh  with  (3=1). 

The  radical  mesh  is  a  special  case  of  a  "graded  mesh"  with  "grading 
function"  g(y)  =  y^,  defined  by 


n  1 1  I  .  , 

"l  '  8  IT  •  ■  '  ' . ”n 

^  n-' 


Remark  2.1:  A  generalization  of  graded  meshes  uses  the  notion  of  density 
(i/>  =  g'  in  the  one-dimensional  case).  It  is  naturally  applicable  to  two 
dimensions.  See  [Hu] . 

2 

We  consider  here  only  a  denumerable  family  of  meshes.  Of  course,  we 
could  consider  a  one  (or  more)  parametric  family  of  meshes  as  well. 


7 


3.  The  geometric  mesh  with  ratio  q,  0  <  q  <  1  (with  respect  to  x  = 
0).  Here, 


(2.7a) 


i  =  l.....m^ 


(2.7b) 


=  0. 


A  more  general  definition  will  be  presented  in  Section  4  for  the  2-d  case. 


, p  )  p  .  2:  1  an  integer,  be  called  the 
^n,n  ^n,j 


degree  vector.  By  V^(3'^,p^)  c  Hq(Q),  we  denote  the  set  of  all  functions  v 
such  that  the  restriction  v|  ^  is  a  polynomial  of  degree  j ’ 


lies  in  T  .  If  p  .  =  p  ,  then  we  speak  about  the  uniform  degree  vector 

p  ^n,j  ^n  ^  - ^ - 

“  ♦  J 

for  which  we  write  p^.  =  V^(3’*',E^)  will  be  called  the  space  of 

elements  of  degree  p  . 

Having  constructed  a  finite-dimensional  subspace  as  above,  we  now  define 
the  finite  element  solution  as  u^  €  satisfying  (2.3)  for  all  v  e  V^. 

Then  we  immediately  obtain,  with  e^  =  u  -  u^. 


(2.8)  e  s  inf  u  -  v  , 

"  n"E  ,,n  "  "E 

v€V 

1/2 

where  =  (B(u,u))  is  the  energy  norm,  which  is  equivalent  to  the 

H^(n)  norm. 

In  this  paper,  we  will  be  most’/  interested  in  the  energy  norm  measure 
of  the  error  in  the  finite  element  solution.  Any  change  of  the  norm  can  (and 
usually  does)  change  the  results  and  conclusions  following  from  them.  In 
practice,  many  different  accuracy  assessments  are  used. 
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2.2  Extension  processes 

In  order  to  increase  the  accuracy  of  the  finite  element  solution,  we 
define  an  extension  process  which  entails  the  construction  of  a  sequence  of 
subspaces  (of  particular  character).  The  classical  extension  process 

uses  where  ,  n  -  1,2,...,  is  a  sequence  of  meshes  with  h^ 

decreasing  (h^  ->0  as  n  -^  «)  and  =  p,  usually  p  =  1,2.  Here  the 
meshes  3"^  could  be  a  sequence  of  uniform,  quasiuniform  or  non-quasiuniform 
meshes.  This  choice  of  is  called  the  h  version  of  the  FEM.CM] 

The  p  version  of  the  FEM  on  the  other  hand  uses  ) ,  where  = 

1  is  a  fixed  mesh  and  the  degrees  p^  ^  oo  as  n  oo.  The  increase  of  p^ 

can  be  uniform  (p_=  p  ).  or  selective  (i.e.,  p  .  are  different  for 
^  n  ^n,j 

different  j).  The  mesh  3"  can  be  essentially  uniform  or  it  may  be  strongly 
refined  in  certain  areas. 

Finally,  the  h-p  version  of  the  FEM  consists  of  V^(3'^,p^)  with  3"^  in 
general  being  refined  and  p^  increasing  with  n.  Obviously,  the  h-p  version 
is  a  generalization  of  the  h  and  p  versions. 

It  is  useful  to  distinguish  between  hierarchic  and  non-hierarchic 
extensions.  Hierarchic  ones  are  such  that  p  V^.  This  assures  that  the 

error  decreases  monotonical ly ,  i.e.,  ll^fi+ll^E  “  ll^n^^E’  spaces  are 

not  hierarchic  then,  in  general,  this  monotone  behavior  does  not  hold.  The 
basis  functions  for  the  p  version  are  in  general  hierarchic. 

In  practice,  of  course,  we  always  use  finite  meshes  with  elements  of 
finite  degree  so  that  for  concrete  computations,  the  distinction  between 
different  extensions  is  blurred.  Nevertheless,  this  distinction  is  important 
not  only  theoretically,  but  also  for  practical  purposes,  since  while  per¬ 
forming  an  actual  computation  it  is  usually  assumed  that  the  results  are  in 
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the  asymptotic  range  of  the  extension  process.  Since,  in  practice,  one  or 
more  steps  of  the  extension  process  are  usually  carried  out,  it  is  here  that 
the  behavior  and  asymptotic  properties  of  various  extensions  play  an  important 
role. 


2.3.  Approximation  theorems  for  quasiuniform  meshes 

The  rate  at  which  the  finite  element  approximations  converge  to  the  true 
solution  will  depend  upon  two  factors:  the  definition  of  each  of  the 
subspaces  (i.e.,  quasiuniform  or  non-quasiuniform  meshes,  choice  of 

etc.)  as  well  as  the  type  of  extension  procedure  (h,  p  or  h-p).  We  now 
present  various  approximation  theorems  which  bound  the  right  hand  side 
of  (2.8),  in  terms  of  the  number  of  degrees  of  freedom  =  dim  V^.  If 
various  methods  are  to  be  compared,  they  should  be  assessed  in  terms  of  the 
total  work  required  to  obtain  a  given  accuracy.  Although  is  one  measure 

of  the  work  required,  it  does  not  give  the  complete  picture,  since  various 
other  factors  like  amount  of  user  time,  special  solution  techniques,  etc. ,  are 
involved  (see  [BE],  [BEM],  [CM]).  With  these  qualifications  in  mind,  we  first 
discuss  the  case  of  the  h-p  version  using  quasiuniform  meshes.  We  will  omit 
the  subscript  n  where  no  confusion  would  occur. 


Theorem  2.1.  Let  =  p^  =  p, 


and 


7 


n 


be  quasiuniform. 


Then  for  k  2:  1 


(2.9)  l(e||  ^  Ch^^~^V"^^"^^l|ui 

H  (n) 

where  p  =  min(k,p+l)  and  C  is  independent  of  n  (i.e.,  of  h,  p)  and  u, 
but  depends  on  k  and  the  quasiuniformity  constant  t. 

From  this,  noting  that  N  =  N(V^)  s:  Ch  ^p,  we  obtain  the  following 
corollaries . 
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Corollary  2.1.  For  the  h  version  (i.e.,  =  p)  we  get 

(2.10)  ||e||  s  Ch^^”^^|u||  i 

H  (Q) 

Corollary  2.2.  For  the  p  version  (i.e.,  =  S')  we  get 

(2.11)  ||e||  i  Cp“^^“^^l|u|l  ^ 

H  (Q) 

The  estimate  (2.9)  is  a  refined  version  of  the  estimate  (2.10)  which 

is  the  one  found  in  many  books  on  the  FEM.  It  shows  explicitly  how  the 

constant  C  in  (2.10)  behaves  with  respect  to  p.  The  two-dimensional 

version  of  this  estimate  (see  also  Theorem  4.1)  was  proven  in  [BSl]  while 

that  of  the  estimate  (2.11)  was  proven  in  (BS2]. 

Notice  that  (2.9)  shows  that  if  the  h  version  is  used  in  a  case  that 

k  <  p  +  1,  increasing  p  further  will  decrease  the  error  by  decreasing  the 
- (k-l ) 

"constant"  Cp  ,  though  it  will  not  affect  the  exponent  of  h.  This 

exponent  of  N  (or  h)  is  often  called  the  asymptotic  order  of  convergence, 
and,  for  example,  for  (2.10)  will  often  be  written  as 

|e|j,  =  OIN-'"-'’). 

The  above  estimates  are  optimal  in  the  following  sense.  For 
any  N,  there  is  a  function  u  e  H  (Q),  possibly  depending  on  N,  such  that 

(2.12)  ||e||  ^  CN'^‘''^^||u|| 

H  (n) 

with  C  >  0  independent  of  u,  N,  the  mesh  J,  and  element  degrees  p. 

(2.12)  follows  from  the  theory  of  n-widths  (see  [Pi]). 

If  k  s  p  +  1,  then  for  the  h  version  there  is  a  u  €  H  (£1)  such  that 
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(2. 13) 


lei  i  C(p)hP(|u((  ^  CN"P 

H  (n) 

with  C  independent  of 

Let  us  remark  that  the  error  bounds  (2.12)  and  (2.13)  must  be 
distinguished  from  asymptotical  expansions  of  the  type 

(2.14)  |e|^  =  C(u)N"^‘^"^^  +  o(n"^‘^"^^ 

which  are  stronger  estimates. 

In  practice,  we  are  usually  interested  in  the  convergence  of  the  FEM  for 
a  specific  u,  namely  the  solution  of  the  problem  of  interest.  For  this 
specific  u,  the  rate  of  convergence  could  be  higher  than  the  rate  for  the 
"worst"  case  (which  is  the  one  addressed  in  (2.9)  and  (2.11)).  For  example, 
function  u  defined  in  (2.2)  belongs  maximally  to  H  (I)  but  the 

convergence  rate  in  (2.9)  can  be  too  pessimistic  since  for  particular 
relations  between  h  and  p,  better  asymptotic  rates  of  the  form  (2.14)  hold 
(see  e.g.  Theorem  2.2).  In  general,  we  must  therefore  distinguish  between  the 
convergence  rate  for  the  worst  case  (which  may  sometimes  only  be  achievable 
for  a  sequence  of  exact  solutions  dependent  on  N)  and  the  convergence  rate 
for  a  particular  case,  especially  when  the  solution  under  consideration  is 
"reasonable".  Consequently,  when  the  accuracy  of  the  FEM  is  experimentally 
studied  by  numerical  determination  of  the  exponent  (k-1)  of  N  in  (2.14), 
it  should  be  kept  in  mind  that  this  only  gives  an  upper  bound  for  the  true 
exponent  in  (2.12),  (2.13)  and  cannot  always  serve  as  an  indicator  of  the 
performance  of  the  method  for  all  cases. 

The  estimate  (2.13)  shows  that  the  maximum  possible  rate  of  convergence 
using  the  h  version  is  always  bounded  by  the  fixed  degree  p  of  polynomials 
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used.  In  contrast,  (2.11)  shows  that  the  p  version  gives  arbitrarily  high 
rates  of  convergence,  depending  oniy  on  the  smoothness  of  the  function  u. 

The  above  estimates  are  based  only  on  the  information  that  u  e  H  (£1). 
If  more  information  is  avaiiable,  then  better  estimates  can  be  obtained.  For 
example,  if  u  is  analytic  on  [0,11,  then  the  error  for  the  p  version 
decays  exponentially. 

Let  us  now  assume  that  the  solution  has  a  singularity  at  a  point,  or 
specifically,  that  u  is  given  by  (2.2).  Then  we  have 


Theorem  2.2.  Let  the  solution  u  be  given  by  (2.2).  Then 

I,  ^  ^  .  min(a-l/2,p-a+l/2)  -2(a-l/2), 

||e||j.  s  C  min[h  ,h  ^  p  ] 

This  gives  immediately 


Corollary  2. 3.  For  the  h  version 
(2.15)  ||e||^  s 


Corollary  2.4.  For  the  p  version 


(2.16) 


i  Cp 


-2(a-l/2) 


=  0(N 


-2 (a- 1/2) 


The  estimate  (2.15)  follows  from  (2.10)  when  Besov  spaces  are 

Z  ,  00 

used  instead  of  fractional  Sobolev  spaces  H  (£2).  It  is  sharp  when 
quasiuniform  meshes  are  used,  in  the  sense  that  there  exists  C  >  0  such  that 


(2. 17) 


Also,  (2.16)  is  sharp  as  well.  Theorem  2.2  is  once  more  a  special  version  of 
the  theorem  for  the  two-dimensional  case  proven  in  [BSl],  [BD]  (where  “  “  ^ 
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has  to  be  replaced  by  a). 

We  see  that  the  p  version  gives  twice  the  rate  of  convergence  as  the  h 
version  (on  quasiuniform  meshes)  when  the  solution  has  an  type  of 
singularity  at  an  endpoint  of  the  element.  If  the  singularity  is  located 
inside  the  element,  then  this  doubling  in  the  rate  does  not  occur.  The  above 
result  about  doubling  of  the  rate  was  first  proven  for  the  two-dimensional 
case  in  [BSK],  where  a  slightly  weaker  result,  showing  the  doubled  rate  of 
convergence  reduced  by  an  amount  0(N  )  (e  >  0  arbitrary),  was  obtained. 

The  dependence  on  c  was  removed  in  [BS2]. 

Let  us  outline  the  reason  for  obtaining  the  rate  of  convergence 
0(N  1/2)  j  Corollary  2.4.  Suppose,  for  simplicity,  we  have  just  one 

element  and  the  domain  is  Q  =  (-1,1)  instead  of  n  =  (0,1),  with  the 
solution  (2.2)  correspondingly  scaled.  Then  if  we  were  to  use  the  standard 
estimate  (2.11)  with  u  being  imbedded  in  the  Besov  space  like  we  did  for  the 
h  version,  we  would  get  an  identical  rate  0(N  1/2)^  (2.15).  It 

turns  out,  however,  that  using  various  orthogonality  properties  of  Jacobi 
polynomials,  one  can  show  that  [GiB] 


(2.18) 


1^  -  %IIe  ^  cp  ‘^lu-  II 

P  vnn) 

k 


where 


(2. 19) 


I  Ir 


(l-t^)‘'(v^’'^^  dt 


so  that  the  usual  Sobolev  norm  in  (2.11)  is  replaced  by  the  weighted  Sobolev 


norm  in  (2.19).  Now  it  may  be  readily  verified  that  although  the  solution 
(2.2)  satisfies  u'  g  H^(n)  for  k  <  a  -  ^  only,  due  to  the  weight  (1-t^)*^ 
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k  1 

in  (2.19),  we  have  instead,  u'  e  for  k  <  2(a  -  •  leading  to  the 

doubled  rate  up  to  an  arbitrary  c  >  0. 

“k 

In  order  to  obtain  a  rate  of  0(p  )  with  k  actually  equal  to 

2(a  -  ^)  without  the  c,  a  more  detailed  analysis  is  required.  Specifically, 
we  expand  u'  in  terms  of  the  Legendre  polynomials 


(2.20) 


u'(x)  = 


I 


a  P  (x) 
n  n 


n=l 


and  choose  u  such  that  u  (±1)  =  0  and  u*  is  the  sum  of  the  first  (p-1) 
P  P  P 

terms  of  (2.20).  Then  we  have,  by  the  orthogonality  property  of 


(2.21) 


u  -  u 


p"E 


2  2 


2n+l 


Here,  a^  may  be  explicitly  computed  by  the  formula 


2n+l 


r  u'P  (t)dt  =  f  (a(l+t)“  ^-2“)  [(1-t^)^]^^ 

J-1  "  2  J-1  2^n! 


dt 


using  the  Rodriguez  form  of  Legendre  polynomials.  Once  again,  the  weight 
(1-t^)^  leads  to  an  "absorption"  of  the  singularity  at  the  end  of  Q. 
Integrating  by  parts  and  using  some  properties  of  the  gamma  function  leads  to 
the  following  estimate,  proved  in  [GiB] 


(2.22) 

with 


1  r-f  1  Cl  C_  (a) 

,  , ,n-l  C(a)  ,,  1  2 

a  =  (-1)  — — -  (1  +  -  +  — - —  +  ••• 

n  2a- 1  n  2 


n 


C(a)  = 


n 


2*  I  r(a)  I  ^sinjr(a-l ) 


(2.21),  (2.22)  lead  to  the  following  theorem,  which  has  been  stated  in  terms 
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of  the  original  domain  =  (0,1)  ([GiBl); 


Theorem  2.3.  For  the  p  version  in  the  case  of  a  single  element,  we  have 


(2.23) 

as  p  00,  where 


^IIe  ^0^“^  2a-l 


(1  +  0(-)) 
P 


CQ(a) 


a(r(a) )  1  sin  nal 


7tv'2a-l 


Let  us  show  a  numerical  example.  Consider  the  cases 
3.5.  In  Table  2.1  we  show  computed  values  of 

Ep  is  the  leading  term  in  the  right  hand  side  of  (2.23). 
results  are  seen  to  be  in  close  agreement  with  (2.23). 


a  =  0 . 7  and 
and 

P  gA 

P 


The  numerical 


a  = 
where 


TABLE  2.1.  The  error  of  the  finite  element  solution. 


a  =  0.7 

a  =  3.5 

p 

E 

r" 

E 

R^ 

P 

P 

P 

P 

1 

4.743E-1 

0 . 9877 

1.021 

0.203 

2 

3.637E-1 

0.0067 

3.402E-1 

4.335 

3 

3.090E-1 

0.9985 

3.093E-2 

4.488 

4 

2.766E-1 

0.9992 

2.379E-3 

1.949 

5 

2.522E-1 

0.9995 

4.760E-4 

1.480 

6 

2.344E-1 

0.9996 

1 . 400E-4 

1.300 

7 

2.204E-1 

0.9997 

5.  154E-5 

1.208 

8 

2.090E-1 

0.9998 

2.210E-5 

1.153 

9 

1.994E-1 

0.9998 

1 . 057E-5 

1.  118 

10 

1.912E-1 

0.9999 

5.495E-6 

1.094 

11 

1 . 840E-1 

0.9999 

3.053E-6 

1.077 

12 

1.777E-1 

0.9999 

1 . 790E-6 

1.064 

13 

1 .722E-1 

1.000 

8.999E-7 

1.054 

14 

1.671E-1 

1.000 

6.978E-7 

1.046 

15 

1.626E-1 

1.000 

4.585E-7 

1.040 

16 


2.4.  Approximation  theorems  for  non-quasiuniform  meshes 

Let  us  now  consider  the  case  when  non-quasiuniform  meshes  are  used  and 
the  solution  is  given  by  (2.2).  We  will  show  that  it  is  possible  to  get  an 
exponential  rate  of  convergence  (in  N)  using  the  h-p  version  of  the  FEM. 
First,  however,  let  us  discuss  the  h  version.  Let  us  show  that  the  best 
possible  rate  is  0(N  ^)  (no  matter  what  mesh  is  used)  i.e.,  the  convergence  is 
always  algebraic.  As  before,  N  is  the  number  of  degrees  of  freedom  and  p 
is  the  degree  of  polynomials  used  in  the  h  version. 

For  simplicity,  we  take  p  =  1.  Let  w(x)  =  u'  (x)  =  ax*^  ^  -  1  on  n  = 


I  .  = 
J 


(Here  m  is  the  number  of  elements).  Since  the  h.’s  in  the  above  are 
n  J 


arbitrary,  we  estimate  the  above  sum  as  follows.  Let 


J 

Then  obviously 


\ 

2 


so  that  since 


(0,1)  we  have 


JeJ 


1 

2 
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Therefore 


J=1  jeJ 


Because  N  =  m  -  1  we  get 
n  ® 


(2.24) 


i®IIe  - 


no  matter  what  kind  of  mesh  is  used.  A  similar  result  holds  for  p  >  1. 

The  following  theorem  shows  that  we  also  have  ^  CN  ^  provided  the 

mesh  is  properly  chosen  (i.e.,  the  above  bound  is  attained). 


Theorem  2. 4.  For  the  h  version,  using  polynomials  of  fixed  degree  (i.e.,  p^  = 
p),  the  radical  mesh  using  m  elements  with 


(2.25) 


0  ~  a-1/2 


is  the  optimal  choice  and  the  following  estimate  holds  with  h  =  - 

m 


(2.26) 


where 


(2.27) 


h-xn  h*' 


C(a,p) 


.  ar(a)  sin  na  r(p-a+l) 
C(a,p)  =  - - - - 

4^  VZ^  r(p+l/2) 


The  notion  of  optimality  in  the  above  theorem  is  understood  in  the 
following  sense.  Consider  the  sequences  of  meshes  defined  by  grading 
functions  g  that  satisfy 

i)  g  €  C°[0,1]  n  C^(0,1) 

ii)  g(0)  =  0,  g(l)  =  1  and  g  is  strictly  increasing. 

Then  the  error  of  the  FEM  with  the  radical  mesh  (g  =  x  )  attains  the 
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minimum  among  ali  grading  meshes  satisfying  i),  ii). 

It  is  better  to  (mildly)  overrefine  the  mesh  than  underrefine  it.  This 
can  be  seen  from  the  detailed  results  proven  in  [GiB] .  Mild  overref inement 
is  advisable  in  two  and  three  dimensions  as  well. 

We  see  that  proper  refinement  leads  to  the  highest  possible  rate  in  the 
h  version,  algebraic  in  N,  which  is  the  same  as  would  be  achieved  if  the 
solution  were  smooth.  We  also  see  that  for  fixed  p  (uniform)  the  best  mesh 
to  use  (for  the  h  version)  is  a  radical  one,  with  the  refinement  strength 
(i.e.,  ^)  being  increased  linearly  with  p. 

Let  us  now  show  that  an  exponential  rate  of  convergence  is  possible 
using  the  h-p  version.  For  simplicity,  we  will  take  a  uniform  degree 
vector.  Let  us  take  the  mesh  to  be  the  geometric  one  defined  by  (2.7)  and  let 


1  (Xl_l  *  X.) 


As  before,  let  w(x)  =  u‘(x),  with  u(x)  given  by  (2.2),  and  let  us  expand 
w(x)  in  terms  of  its  Taylor  expansion  about  the  point  x^  so  that 

( i )  k  ^  f  i )  k 

(2.28)  w(x)  =  w.  (x)  +  r.  (x)  =  )  a'  ^(x-x.  )^  +  V  a,  ^(x-x.  )^ 

i,p  i,p  L  k  1-1/2  L  k  1-1/2 

k=0  k=p 

where  w.  (x)  is  the  Taylor  polynomial  of  degree  p  -  1  and  r.  (x)  is 

i.p  j  ^  j  or- 

the  remainder.  Using  (2.7a)  (2.7b)  we  see  that 


x.+x.  ,  .  , 

1  1-1  m-i.l+q. 

’‘i-.XE  =  —2—  =  " 


x.-x.  ,  .  , 

1  1-1  m-i,l-q, 

- 2 -  "  ^  ^  2  ^ 


Now  for  i  =  2,3,...,m,  the  function  x**  ^  is  analytic  in  the  closed  disk 
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for  i  i  3) 


with  center  x,  .  _  and  radius  q^'"  -  q^)  (=  x. 

(m-i ) (a-1 ) 


_  _  q  -  q  ;  1=  X.  “  X.  „ 

1-1/2  2  1-1/2  1-2 


and  is  bounded  by  Cq 


there.  Hence,  applying  the  Cauchy  theorem, 


where  |p|  =  |2q+l|  ^  <  1,  and  C  is  independent  of  m,  i,  k.  Hence 


(2.29) 


,  ,  ,  ,  w  ^  r-  (m-i)  (a-1)  p 

I w(x)  -  w^  p(x) I  i  Cq  pT . 


Using  this. 


(2.30)  =  inf 


J, 


r  -  w  dx  s 


re^  ,  •’  I . 

p-1  1 


1  w  .  -  w  I 

Jt  '  1>P 


dx 


'I. 

1 


2p  (m-i ) (2(a-l )+l ) . 
Cp  q 


Therefore 

(2.31) 


m 


i=2 


Next,  we  note  that 


(2.32)  inf  r  |r-w|^dx  s  f  |w{x)|^dx  =  Cq 

Jt  Jt 


re^’  ,  •'I, 

p-1  1 


(m-1 ) (2a-l )  m 

-  ‘=‘>1 


for  some  0  <  p^  <  1.  Choosing  p  =  xm  and  noting  that  N  =  m(p+l)  we  see 
using  (2.31),  (2.32)  that 


i=l 


-yN 


1/2 


for  some  y  >  0  and  C  depending  on  q,  a. 


Now  let  r,  denote  the  minimizer  of  E,  in  (2.30).  Then  obviously 
i  1 
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0. 


(r.-w)dx  = 
1 


Hence  defining  R(x),  0  <  x  <  1,  R(x)|j  =  r^(x)  we  can  define 


<j(x)  =  [ 

Jn 


R(t)dt 


and  get  wClJ  =0.  Hence  o)  €  V  (ST.p)  such  that 


||u-a)||  s  Ce 

h"(I) 


-yN 


1/2 


and  the  exponential  rate  is  proven.  See  [GiB]  for  a  more  precise  analysis. 

We  have  used  p  uniform.  We  could  also  analyze  the  case  for  p 
non-uniform.  In  fact,  the  optimal  rate  is  obtained  by  choosing  a  non-uniform 
degree  vector  as  described  in  the  following  theorems. 


Theorem  2.5.  Let  [a]  denote  the  integral  part  of  a.  Then  the  optimal  mesh 
for  the  h-p  version  is  the  geometric  one  with  ratio  =  (v^  -  1)^  =  e  ^ 

»  0.17  and  linear  degree  vector  p  =  ^PjK  Pj  =  +1  with  s^  =  2a  - 

1.  With  this  choice. 


(2.34) 


I  II  ^  \  v'(a'l/2)N  ,  -1.7626V'(a-l/2)N 

e||„  s  C(a)q  =  C(a)e 


Let  us  be  more  precise  about  the  notion  of  optimality.  We  have  the 
following  theorem. 


Theorem  2.6.  For  any  mesh  H ,  and  any  degree  vector  p,  we  have 


lellp  -  C(a) 


✓(a-l/2)N 

^0 _ 
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Comparing  the  lower  bound  given  in  Theorem  2.6  and  the  upper  one  in  Theorem 
2.5.  we  see  that  the  geometrical  mesh  with  linear  degree  vector  gives  the 
maximal  possible  exponential  rate.  This  is  the  justification  for  calling  the 
mesh  optimal. 

For  any  given  geometrical  mesh  with  ratio  q  there  is  an  optimal  p, 
which  is  linear,  with  slope 


(2.35) 


Sfllq) 


log  q 

log  r  ’ 


r 


1-v/q 
^  ■ 


For  this  choice,  the  error  is 


(2.36) 


s  C(a.q)e”'^^“"^''^^^  vTTog"VT^g~r 


This  shows  that  for  small  q  (i.e.  strong  refinement)  the  optimal 
degree  vector  is  less  uniform. 

In  Theorems  2.5  and  2.6,  we  have  addressed  the  case  when  the  degree 
vector  was  non-uniform.  For  uniform  degree  vector,  the  situation  is  slightly 
different. 


Theorem  2.7.  Let  q  be  the  ratio  of  the  geometric  mesh  with  m  elements  and 
p  =  s^m  where  s^  is  as  in  (2.35).  Then 


where 


C(a.q)e‘'^^“'^^2’^  ^-<r/2 


(T  =  min(2a-l,a). 


Substituting 

(2.37) 


the  optimal 


2 

q  =  q^  =  -  1)  in  the  above, 

,  -1.2464v'(a-l/2)N  ,,-o-/2 
£  C(a)e  N 


we  get 


Hence  we  see  that  the  exponential  rate  is  smaller  in  this  case. 
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In  practice  (in  higher  dimensions)  one  often  uses  a  geometric  mesh  with 
a  fixed  number  of  elements  m  and  determines  p  more  or  less  adaptively. 

For  every  m,  the  convergence  will  have  two  parts.  In  the  pre-asymptotic 
portion,  the  rate  is  exponential,  while  in  the  asymptotic  region,  one  obtains 
the  algebraic  rate  of  the  p  version.  If  the  result  of  a  computation  lies  in 
the  algebraic  part,  then  this  Indicates  the  mesh  was  not  sufficiently  refined. 

If  a  uniform  degree  vector  p  is  to  be  used,  then  it  turns  out  that  the 
geometric  mesh  does  not  give  the  best  rate  of  convergence,  i.e.  (2.37)  is  not 
optimal.  We  have  seen  that  for  fixed  uniform  degree  vector  (i.e..  -  h 

version),  the  radical  (and  not  the  geometric)  mesh  is  'optimal  .  If  we 
consider  a  sequence  of  uniform  degree  vectors  and  perform  the  h  version  for 
each  of  them  with  the  corresponding  optimal  radical  mesh,  then  the  envelope  of 
the  resulting  error  curves  leads  to  an  h-p  version  which  gives  better  results 
than  the  one  in  Theorem  2.7.  We  have  the  following  theorem. 


Theorem  2.8.  There  is  an  h-p  version  with  uniform  degree  vector  p  which 
has  the  following  estimate 


(2.38) 


C(a)  -(4/e)V(a-l/2)N 

E  "  „(a-l/2)/2  ® 

N 


C(a)  -1.4715v'(a-l/2)N 

,,(a-l/2)/2  ® 

N 


As  N  -a  00,  the  meshes  (inside  (0,1))  tend  to  the  geometric  mesh  with 
-4/e^ 

ratio  q  =  e  ~  0.5820  and  the  relation  between  the  degree  p  and  the 

number  of  elements  m  tends  to  a  linear  one. 


P 


Various  numerical  example  illustrating  these  results  may  be  found  in 
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[GiB] . 


2.5.  The  rate  of  convergence  in  the  norm 

So  far,  we  have  discussed  the  h,  p  and  h-p  versions  in  terms  of  their 
convergence  in  the  energy  norm.  An  important  difference  between  the 
h  version  on  the  one  hand  and  the  p  (and  h-p)  version  on  the  other  that  is 
apparent  from  the  above  is  that  the  p  (and  h-p)  version  leads  to  different 
asymptotic  rates  of  convergence,  depending  upon  the  nature  of  the  solution 
(for  example,  the  "doubling"  effect  for  singularities).  In  contrast, 

the  only  way  that  the  h  version  could  lead  to  a  rate  better  than  O(h^)  for 
a  function  f(x)  is  if  f  were  a  polynomial  of  degree  ^  p.  This 
difference  leads  to  an  interesting  effect  in  the  error  in  the  and  other 

lower-order  norms. 

Let  us  discuss  the  error  in  the  norm  for  the  p  version.  Suppose 

we  first  consider  the  case  where  u  €  H  (Q).  Then,  by  the  usual  duality 
argument ,  we  may  show 


(2.39) 


SO  that  using  Corollary  2.2, 


I®IIl  (£2)  -  Ill'll  k 

2^  ^  H  (£2) 


In  fact,  the  usual  duality  gives  the  following  theorem. 


Theorem  2.9.  Let  =  p^^  =  p  and  be  quasiuniform.  Then  for  k  i  1 , 

(2.40)  ||e||  ^  ChV'^llul 

2^  ^  H  (£2) 
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where  m  =  min(k,p+l). 

This  result  is  optimal  in  the  sense  that 


u€H^(n) 

M  k 

H^(n) 


veV  2 


If  we  now  consider  the  special  singular  solution  of  the  form  (2.2) 

(scaled  on  to  n  =  (-1,1)),  then  Theorem  2.9  predicts  a  rate  of 

min(a+l/2,p+l ) ,  u  w  i.  •  , 

0(h  )  for  the  h  version  which  is  the  best  error  possible 

with  a  quasiuniform  mesh,  and  so  is  again  optimal.  For  the  p  version  (using 

a  single  element),  we  may  now  use  (2.39)  with  (2.23)  to  obtain  an  0(p  ) 

rate.  The  optimal  approximation  estimate  is  essentially  obtained  by 

taking  a  =  a  +  1  in  Theorem  2.3,  which  gives 


(2.41) 


-2(a^|) 

inf  |u-v||j^  =  Cp  (l+0(p  )) 


Hence  the  estimate  through  duality  does  not  give  the  optimal  rate  (2.41). 

-Z«-1 

An  improved  estimate  of  0(p  )  may  be  obtained  using  the  fact  that 


(2.42) 


P  (t)dt  =  (P  ^,(x)  -  P  ,(x)) 

p  2p+l  p+1  p-1 


Then  by  (2.20),  the  definition  of  u  ,  and  (2.42), 

P 


(2.43)  u  -  u 


<»  a  pc 

I 


(t)dt 


2p+l  p-1 


p  ^  y  — 

2p+3  p  L  l2r+l 


2r+5  r+r 


25 


Def ine 


b 

r 


a 


r 

2r+T 


r+2 

2r+5 


Then,  using  (2.22),  we  obtain 


(2.44) 

Now  by  (2.43), 

(2.45)  II 


C,(a)  C2{a) 


^  ,  ,r-l  C(a)  (,  ^2'  '  1 

r  V  r  ' 


u  -  u 


1_  +  2 


r  1-2  2 

+1  ^  Z  2r+ 


(2p*l)2  2P-1  ,2p,3,2  2p 

Using  (2.44),  we  that 


r=D 


I 

r=p 


C(a) 

4a+2 

P 


(1 


+  0(p"^)). 


Also,  (2.22)  may  be  used  to  estimate  the  first  two  terms  in  (2.45),  yielding 


(2.46) 


u  -  u 


'L2(^) 


C(a) 

2a+l/2 

P 


(1 


+  0(p"^)). 


This  estimate,  which  is  sharp,  improves  upon  the  rate  through  duality,  but 

still  does  not  give  the  optimal  rate  of  (2.41). 

Hence  we  observe  the  following  interesting  difference  when  the 

projection  of  the  function  (2.2)  is  taken  in  the  energy  norm  (equivalent  to 

the  H^(n)  norm).  For  the  h  version,  the  convergence  in  the  ^2^^^  norm  is 
min(a+-, p+1 ) 

0(h  )  which  is  optimal.  For  the  p  version,  this  convergence 

-2ol-- 

2 

rate  is  higher,  0(p  ).  However,  the  best  approximation  is 

1 

0(p  ),  which  is  even  higher.  As  a  result,  the  energy  norm  projection 

does  not  give  the  optimal  rate  in  the  L^ID)  norm  for  (2.2)  (see  [Er]).  Let 
us  mention  that  in  two  dimensions,  for  non-convex  domains,  the  h  version  also 
fails  to  give  the  optimal  rate  (see  Section  4.2). 
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2.6.  The  pollution  problem 


Let  us  now  discuss  another  difference  between  the  h  and  p  versions,  the 
phenomenon  of  pollution.  This  is  said  to  occur  when  unsmoothness  of  the 
solution  u  in  one  area  influences  the  accuracy  of  u^  in  a  different  area, 
i.e.  the  accuracy  of  u^  is  not  governed  by  the  local  behavior  of  the 
solution  u.  Consider  once  again  the  approximation  of  u  given  by  (2.2) 
(scaled  to  (-1,1)).  Then  u  is  singular  at  x  =  -1  but  analytic  everywhere 
else.  Suppose  u^  is  the  solution  using  the  h  version  (say  with  p  =  1). 

Then  it  is  easily  proven  that  u^  is  simply  the  linear  interpolant  of  u  (an 
analogous  result  holds  for  higher  p).  This  obviously  implies  that  the  error 
|(u-u^)(0)|,  for  instance,  will  only  depend  upon  the  behavior  of  u  over  some 
I-  =  (-5,5),  5  <  1.  In  fact,  this  error  will  converge  asymptotically  at 

O 

2 

the  optimal  rate  of  0(h  )  for  any  a  >  1/2.  Also,  for  fixed  6. 

||u  -  u  I  for  r  =  0,1  will  similarly  converge  at  the  optimal  rate 

H^'d^) 

0(h^  ^),  even  if  1/2  <  a  <  3/2.  We  say,  therefore,  that  the  h  version  is 

free  from  pollution.  The  h-p  version  also  turns  out  to  be  free  from 

pollution,  so  that  |  (u-u  )(0)|,  ||u  -  u  ||.  .,  ||u  -  u  ||  all  converge 

n  n  n  ) 

O 

at  the  optimal  exponential  asymptotic  rate  determined  only  by  the  analytic 

behavior  of  u  over  I  , ,  5'  close  to  5,  and  not  by  the  singular  behavior 

o 

at  X  =  -1. 


The  p  version,  however,  will  exhibit  pollution.  Suppose  we  use  a  single 
element,  then  with  u  s  u  ,  ||u  -  u  ||  ,  r  =  0,1,  will  be  no  better  than 

"  P  '  Ph^(I^) 


the  rates  given  in  Sections  2.3,  2.5  for  llu  -  u 


D "  r 
P  H  (I) 


,  even  though  u  is 
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analytic  over  I  .  More  precisely,  let  us  estimate  | (u-u  ){0)|.  Noting  that 


,m 


,o>  =  (-1,”  .  iiii- 1. 

2m  ( 2m ) ! ! 


Vn  Vin 


1  +  _  +  •  . 
m 


we  get  for  p  =  2k+l  (using  (2.22),  (2.43)), 

C  -1 


(u-Up)(01  = 


(1  +  0(p  ))  +  J] 


oo  (-1 )  b 


m=k 


Vn  (m+1) 


2m+l  -1, 

(1  +  0(m  ) 


1/2 


Using  (2.44)  and  the  fact  that  the  terms  in  the  sum  are  of  alternating  sign, 
this  gives 

-  TSTI72  »  * 

P 

The  above  estimate  clearly  shows  the  pollution  effect,  since  in  the  absence  of 
pollution,  the  rate  should  be  exponential  because  u  is  analytic  over  . 
The  reason  this  pollution  occurs  should  be  clear:  unlike  the  h  and  h-p 
versions,  there  is  no  "buffer"  of  elements  now  to  isolate  the  singularity  from 
the  rest  of  the  domain.  For  more  results  on  the  pollution  problem  in  1-d,  in 
the  context  of  the  h-version  as  well,  see  [W2] . 


3.  THE  MODEL  PROBLEM  AND  ITS  REGULARITY 

In  Section  2  we  have  analyzed  a  simple  one-dimensional  problem.  In  this 
section  we  will  address  the  basic  properties  of  the  analogous  problem  in  two 
and  three  dimensions.  We  will  concentrate  on  model  problems  which  are  typical 
in  the  field  of  computational  (structural)  mechanics.  Mathematically,  these 
problems  are  described  by  elliptic  partial  differential  equations  with 
piecewise  analytic  input  data.  Hence  this  class  of  practical  problems  is 
relatively  narrow  and  well-defined,  leading  to  several  important  common 
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properties  which  can  be  expioited  in  their  numerical  treatment. 


3.1.  The  model  problem 

2 

Let  us  consider  £2  c  R  to  be  a  curvilinear  polygon  with  piecewise 

M  _ 

analytic  boundary  an  =  U  ,  where  P^^  is  an  open  arc  and  P^^  is  analytic 

i=l 

with  endpoints  . ^\l+l  ~  will  denote  the  vertices 

with  internal  angles  w . w.,  (0  <  cj.  s  2n) .  (If  u.  =  2n  then  we  have  a 

slit  domain,  where  the  boundary  is  taken  to  be  two-sided. )  an  will  consist 

of  two  parts,  P^  =  U  P^  being  the  Dirichlet  boundary  and  P^^  =  an  -  P^^ 

ieZ) 

the  Neumann  boundary.  It  should  be  noted  that  the  boundary  could  be  smooth 
everywhere,  with  =  rr  and  the  vertices  corresponding  to  points  in 

Pq  n  Pi^.  We  restrict  ourselves  to  simply  connected  domains  in  order  to 
simplify  the  exposition.  An  example  of  such  a  domain  and  the  notation  used  is 
shown  in  Figure  3.1. 

3 

In  three  dimensions,  n  c  R  will  be  a  bounded  domain  with  piecewise 

analytic  boundary  dQ  consisting  of  faces  P. ,  i  =  1,2,...,M  which  are 

3  3 

curved  polygons  in  R  ,  Joined  by  edges  ,  i  =  1 ,  .  .  .  ,  n^  (curves  in  R  ) 

and  vertices  A.,  i  =  1 . n  .  P^,  and  P.,  will  now  be  the  union  of  some 

1  w  D  N 

faces. 

Analogously  to  (2.1)  we  will  now  consider  a  second  order  elliptic  system 

of  m  differential  equations  for  the  unknown  vector  u  =  (u, . u  )  which, 

1  m 

when  cast  into  the  variational  form  B(u,v)  (analogous  to  (2.3)),  is  given  by 


(3.1) 


B(u,v) 


.  /■  m  t  9u  .  5v.  .. 

=  (  I  I 


i,J=l  k,f=l 


where  t  =  2  or  3  and  ~  ^Jifk  analytical  functions  on  Q 
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satisfying  an  ellipticity  condition  that  ensures 


(3.2) 


Let 


B(u,v)  >  yju|^, 

H  (Q) 


y  >  0. 


hJcQ)  =  {v  e  H^(Q),  V  =  0  on  r^^}. 


A,, 


Further,  let 


F  be  a  continuous  functional,  defined  on 


of  the  form 


(3.3) 


F(v) 


=  f  f^*vdx  -t-  r 


ii- 


vds. 


where  t  is  analytic  on  n  and  fi  is  analytic  on  every  F^.  c 

Moreover,  let  g  e  C^(F_)  and  be  analytic  on  each  f.  c  F,^.  Then  our  model 
®  D  1  D 

problem  consists  of  finding  u  €  H^(n)  satisfying 


(3.4) 


B(u,v)  =  F(v),  V  V  e  Hp(n) 


30 


(3.5) 


u  =  g  on 

For  Fp  *  0  condition  (3.2)  ensures  a  unique  solution  for  (3.4),  (3.5).  For 
Fp  =  0,  we  need  an  additional  solvability  condition  on  F  and  n. 

We  define  the  energy  norm  |u|^  =  (B(u,u) )  ^'^^,  equivalent  to  the  H^(Q) 
norm.  By  restricting  our  formulation  to  the  above  class  of  problems  with 
piecewise  analytic  input  data  of  the  type  assumed,  we  ensure  that  the  solution 
will  possess  certain  special  features. 

Let  us  now  describe  two  concrete  examples  of  (3.4),  (3.5)  in  which  we 
shall  be  interested.  The  first  problem  is  a  scalar  one  (m  =  1)  on  12  c  r'' 
with 


(3.6) 


B(u.v) 


'  I  ( I  ‘'ilr  -fe) 

“  i=l  ^  ^ 


>  0. 


For  s  1,  we  get  the  Poisson  equation.  Most  of  our  cited  results  are  in 
two  dimensions  but  the  three-dimensional  problem  will  be  mentioned  as  well. 

In  Section  5  we  will  be  interested  in  the  2-d  case  that  k^  are  very 
different  in  magnitude.  Obviously,  in  this  case  the  bilinear  form  becomes 
degenerate  as  ->  0.  Such  degeneracies  occur  in  various  practical 

problems  and  their  numerical  treatment  will  be  discussed  in  the  context  of  the 
robustness  of  the  method  (the  locking  problem). 

Our  second  example  is  the  typical  model  problem  of  two-dimensional 
elasticity.  We  will  restrict  ourselves  here  to  the  case  of  isotropic 
homogeneous  materials  (plane  strain)  for  which  m  =  2,  t  =  2  and 


(3.7; 


I„  [  I 


IJ 


.(u) 


IJ 


[v) 


1-21’ 


(div  u  div 


v) 


i.  j  =  l 
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is  the  Poisson 


where  E  >  0  is  the  Young’s  modulus  of  elasticity,  0  ^  ^ 

ratio  and 


(3.8] 


/  \  1 
e. . (u)  =  ^ 

ij  2 


rdu , 
3x 


3u 

-  ^  • 
3x. 

J  i-* 


For  ^  ^  2'  bilinear  form  becomes  degenerate  as  in  the  first  example. 

Nevertheless,  there  is  a  large  difference  in  these  two  cases.  In  the  first 
one,  the  solution  becomes  less  smooth  as  the  form  degenerates  and  in  the 
limiting  case  the  problem  is  not  elliptic.  In  the  second  example,  the 
degeneracy  does  not  influence  the  smoothness  of  the  solution  and  the  limiting 
case  is  still  an  elliptic  problem.  See  e.g.  [AFl]  for  a  general  analysis  of 
the  degeneracy  with  the  smoothness  being  preserved. 

In  three-dimensional  elasticity  the  formulation  is  analogous  to  (3.7), 

(3.8). 


3.2.  Regularity  of  the  two-dimensional  problem 


for  Q  c  R 
=  an.  Then 


Let  us  now  analyze  the  regularity  of  the  solution  of  our  model  problem 

.2 


First  assume  that  dQ  is  smooth  and  that  either  F,,  =  5f2  or  F^ 

N  U 


(3.9) 


lui 


H^(n) 


-  k-2 

H  (n) 


l|gi  k-l/2  ^ 

H  (Fp) 


If  the  data  (i.e.,  Sfl,  g,  Fi)  are  not  sufficiently  smooth  (at  vertices) 

or  both  parts  F|^  and  F^^  are  present,  then  (3.9)  does  not  hold.  (Note  that 

the  analyticity  of  g,  Ii  is  only  on  F. . )  However,  under  our  assumptions, 

M  ^ 

will  be  analytic  on  n  -  U  A.  and  in  the  neighborhood  of  a  vertex  A.  ,  for  the 

i  =  l  ^  ^ 

scalar  equation  (m  =  1),  u  will  have  the  form  (see  [CS],  [Gr2],  [Gr4] ,  [Da]) 
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C4. 


sin  OL.ip 
J 

cos  a.w 
J 


for  the  Dirichlet  or  mixed  problem 
for  the  Neumann  problem. 


The  mixed  problem  here  means  Dirichlet  conditions  for  <p  =  0  and  Neumann 
conditions  for  (p  =  o>. 

The  number  T  in  (3.10)  is  different  from  0  in  the  following  cases: 
i)  The  coefficients  of  the  differential  operator  are  not  constant; 
ii)  the  operator  contains  lower  order  terms; 
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iii)  the  arcs  F.  are  curvilinear. 

1 

In  the  case  of  the  Laplacian,  the  singular  term  due  to  the  curvature  of 
(with  J  =  1,  t  =  1)  may  be  less  regular  than  the  second  term 
(j  =  2,  t  =  0).  This  happens  if 

+  1  “2  ^  ^  " 


for  all  three  boundary  conditions. 

The  coefficients  c.  have  physical  meaning  and  often  are  called 

J  S  L 

stress  intensity  factors.  For  details,  see  [CS] .  We  see  here  that  the 
solution  is  more  complicated,  but  the  singularity  of  the  solution  is  analogous 
to  (2.2),  the  singularity  in  the  one-dimensional  setting. 


Taking  the  correct  number  of  terms  in  (3.10),  we  can  write 


where  is  a  smooth  cut-off  function  in  the  neighborhood  of  A^.  Then  we 

have 


(3. 12) 


=  c(|| 


I  jf-2 

h*'  ‘"(n) 


|g|l  V_1 


1  N 


nk  3/2^j,  j 

1 


i.e.,  the  smoothness  of  u^  is  analogous  to  that  when  the  domain  and  data  are 
smooth  as  in  (3.9).  Note  that  although  the  exponents  in  (3.10)  depend 

continuously  on  the  angles  the  value  of  S  does  not.  As  a  result,  the 

constant  C  in  (3.12)  can  be  very  large  for  angles  close  to  the 

exceptional  values  for  which  S  =  1  instead  of  0.  This  problem  can  be 
avoided,  however,  by  adjusting  the  expansion  of  u. 
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So  far,  we  only  discussed  the  scalar  equation,  which  was  represented  by 
the  example  of  Laplace’s  equation.  The  case  of  a  system  (like  the  elasticity 
problem)  is  more  complicated  but  very  analogous.  See,  for  example,  [Sa] .  The 
decomposition  (3.10)  still  holds  with  being  a  vector.  In  contrast 

to  the  previous  case,  may  now  be  complex  (and  there  are  more  exceptional 

cases  than  before). 


3.3.  Countably  normed  spaces 

The  regularity  of  the  solution  in  our  setting  of  piecewise  analytic 

input  data  can  be  very  advantageously  described  in  the  form  of  countably 

normed  spaces.  (For  details,  see  [BG2] , [BG3] . )  For  5  <  1,  let  3.=  {x  s  Q, 

o 

r  <  5>  and  let  $  =  r^,  for  0  <  /3  <  1 .  By  Ho’^(S-),  m  ^  t  >  0 

p  P  o 

integers,  we  denote  the  completion  of  the  set  of  all  infinitely  differentiable 
functions  under  the  norm 


“/3  ^^5' 


H  (S^)  2  6 


where 


ID^I^ 


I 


a  ,2 


a  =m 


and 


H^’°(s  ; 

P  o 


(U 

I 


k=0 


1‘^/3+k!^  (S  )■ 

^  o 


I  t  I 

Further,  let  us  introduce  $-(S.),  1^0  consisting  of  those  u  6  H,,’  (S^. 

p  o  13  S 

for  which 


Z  o 
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V.  =  t,l+\ .  C^l,  dil  independent  of  k. 

2 

Finally,  introduce  the  space 

p  o 

e^cs^)  =  {u6h2-2(s^)  I  iD'^ujix)  i 

We  then  have 

$q(S  )  c  io^(S-)  c  (S  ),  e  >  0  arbitrary. 

3  3  3  o  3+e  5 

Under  the  assumptions  made  above,  we  have  the  following  theorem. 


Theorem  3. 1  [BG2]  [BG3] .  Let  u  be  the  solution  corresponding  to  a  second 
order  (scalar)  operator. 

1)  Let  n  be  a  polygon  (i.e.,  F.  are  straight  lines).  Then  the 

m 

solution  u  is  analytic  in  Q  -  U  A. ,  and  in  the  neighborhood  of  any  A. , 

i=l  ^  ^ 

2 

we  have  u  e  B^(S^),  V3  <  5^  where  0  <  3  <  1  depends  on  the  angle  . 

ii)  For  curvilinear  polygons,  u  is  analytic  on  Q  -  UA^  and 

u  €  eo^S.),  0  <  3  <  1- 

3  o 

The  characterization  given  in  Theorem  3.1  is  very  useful  in  practice 
because  it  allows  us  to  construct  the  h-p  version  so  that  its  convergence  rate 
is  exponential,  without  knowing  the  decomposition  form  (3.10)  exactly.  For 
analogous  results  for  problem  (3.7)  we  refer  to  [GoB]. 


3.4.  Regularity  of  the  three-dimensional  problem 

In  three  dimensions,  the  situation  is  more  complicated.  The  solution  is 
again  analytic  on  n\(edges  and  vertices).  There  is  a  singular  behavior  in 
the  neighborhood  of  the  edges  (called  edge  singularities)  where  the  solut.ion 
is  smooth  along  the  edge  and  singular  in  the  direction  perpendicular  to  the 
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edge.  In  addition,  near  the  vertices,  the  solution  has  an  essentially  radial 
vertex  singularity,  as  well  as  transitional  edge-vertex  singularities  that 
appear  in  cones  around  the  edges.  The  stress  intensity  factors,  i.e  analogues 
to  c  .  ,  are  now  functions  along  the  edges  which  are  smooth  but  have  singular 

J  s  t 

behavior  in  the  neighborhood  of  the  endpoints  (the  vertices).  These  stress 
intensity  factors  may  be  bounded  or  unbounded  in  dependence  on  the  geometry. 

In  addition,  there  are  stress  intensity  factors  for  the  vertex  singularity. 

As  an  illustration,  let  us  present  some  regularity  results  from 
[Grl]  [Gr2l,  (which  were  used  in  [D2l  for  the  following  model  problem;  Given 
f  e  c'”(fl),  find  u  such  that 


(3. 13) 


-Au  =  f  in  Q 


u  =  0  on  r. 


where  12  is  a  polyhedral  domain  in  IR  . 

First  we  describe  the  singularity  close  to  n ,  an  open  segment  of  an 
edge  Tq  of  Q  assumed  to  be  along  the  x^-axis.  Let  the  dihedral  angle  at 
be  cj,  then  this  singularity  is  characterized  only  by  the  values  = 
jjt/u).  Then,  in  a  neighborhood  G  of  y  (such  that  G  does  not  intersect 
any  corners  or  edges  other  than  y^)  we  have,  analogously  to  (3.10)  (for  f  e 
C“(Q) ) , 


(3. 14) 


III  '^jst^^G^’^jst^^^"  " 

j=l  s=0  t=0 


where  (r.0,x_)  are  local  cylindrical  coordinates,  c  .  €  C  (y)  and  i/i . 

3  jst  ^jst 

are  analytic  in  0.  As  before,  S  =  0  for  the  problem  (3.13)  except  for 


integral  values  of  a^.  for  which  S  =  1. 
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Next,  we  consider  the  behavior  of  u  in  the  neighborhood  of  a  vertex. 
Let  Aq,  the  vertex,  be  situated  at  the  origin  and  be  an  endpoint  of  the  edge 
Tq.  Let  K  be  the  intersection  of  Q  and  a  solid  circular  cone  with  vertex 
Aq  and  axis  along  which  does  not  intersect  any  edge  having  an  end  point 

at  Aq  other  than  Let  =  K  a  where  is  a  small  ball  of  radius 

Pq  about  the  origin.  Also,  let  K  be  the  infinitely  extended  cone  coinciding 
with  n  in  a  neighborhood  of  the  origin  and  let  R  =  K  a  S,  where  S  is  the 

unit  sphere.  Let  0  <  ^  be  the  eigenvalues  of  the 

Laplace-Beltrami  operator  on  R  with  homogeneous  Dirichlet  boundary  data  on 

/  r  ~ 

5R,  let  ttj  be  the  coefficients  (corresponding  to  the 

singular  coefficients  above)  for  this  operator  for  the  singularity  in  a 

00 

neighborhood  of  a  vertex  with  angle  w.  Then,  with  f  e  C  (Q).  we  have  on 

'^0 

=  I  I  I  (I=]sa0 

j=l  s=0  t=0  i=0 

-1/2+^. 

diP  ^  +  Uq, 

i=0 

where  (p,e,^)  are  the  local  spherical  coordinates,  with  origin  at  A^,  and 
with  0=0  corresponding  to  edge  The  functions  ^Jst’  ®i 

Uq  are  all  smooth,  their  exact  regularity  being  stated  in  Lemma  4.3  of  [D2]. 
Here,  the  first  term  is  a  transitional  edge-vertex  singularity  corresponding 
to  edge  y^  and  vertex  A^,  while  the  second  term  is  the  vertex  singularity  of 
radial  type  which  will  be  present  throughout  R. 

EquationsO.  14) ,  (3.15)  then  characterize  the  two  types  of  singular 


-1/2+p.  a.+t 
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behavior  (near  edges  and  vertices)  for  three-dimensional  problems.  For  the 
singular  behavior  of  the  elasticity  problem  we  refer  to  [Gr3],  [Pe]. 

We  see  that  in  two  dimensions  (for  the  Laplace  equation),  the 
singularity  depends  only  on  Wj,  the  interior  angles  of  the  domain  and  upon 
the  boundary  conditions.  In  three  dimensions,  the  situation  is  similar  for 
edge  singularities,  but  not  for  the  singularities  near  vertices,  which  depend 
upon  the  Laplace-Beltrami  eigenvalues.  This  makes  the  a  priori  determination 
of  the  exact  form  of  the  singularities  more  difficult.  In  general,  for  any 

vertex  ,  i  =  1,2 . n^,  we  now  define  to  be  the  portion  of  the 

unit  sphere  subtended  by  the  infinite  cone  which  coincides  with  n  in  a 
neighborhood  of  A^.  Then,  in  addition  to  the  dihedral  angles  w^,  J  = 
l,2,...,n^,  the  regularity  also  depends  upon  the  coefficients  ^  ® 
where  ^  ^  2  ~  eigenvalues  of  the  Laplace-Beltrami 

operator  on  with  homogeneous  boundary  conditions  on  3S^.  These 

eigenvalues  are  directly  related  to  the  Steklov  eigenvalue  problem  on  K^. 

See  [BPA] . 

Let  us  mention  that  we  can  also  describe  the  regularity  of  the 
solution  using  weighted  spaces,  as  was  done  in  Theorem  3.1  for  the  two- 
dimensional  case.  Unlike  the  isotropic  weighted  spaces  used  in  the 
two-dimensional  case,  however,  the  natural  spaces  in  three  dimensions  to 
describe  the  different  types  of  singular  behavior  in  the  neighborhood  of  edges 
and  vertices  are  anisotropic .  We  will  not  elaborate  on  this  here.  For 
regularity  and  error  estimates  in  terms  of  countably  normed  spaces  in  3-d,  we 
refer  to  [BG4]. 
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4.  FINITE  ELEMENT  SPACES  AND  BASIC  APPROXIMATION  RESULTS 

We  will  proceed  here  analogously  to  the  one-dimensional  case.  Assume 

for  now  that  g  =  0.  Then  as  in  Sections  2  and  3  we  will  consider  c 

H^(Q)  to  be  the  finite  element  space  of  dimension  (number  of  degrees  of 

freedom).  The  finite  element  solution  u  e  is  then  such  that 

n 

B(u  ,v)  =  F(v)  V  V  €  V^. 

n 

4.1.  The  finite  element  spaces 

We  will  now  restrict  ourselves  mainly  to  the  two-dimensional  problem  and 
define  the  meshes  and  elements  analogously  as  before. 

Let  S  =  (-1,1)^  and  T  =  {C.t)|0  <  n  <  (C+Dv'a  ,  -1<C£0,  0  <  t)  < 

(l-^)v'3  ,  0  s  ^  <  1}  be  the  standard  square  and  standard  triangle 

respectively.  The  mesh  3"^  is  now  a  partition  of  Q  into  the  open 

curvilinear  triangles  or  quadrilaterals  denoted  by  J  ~  . % 

such  that  to  every  associate  an  invertible  mapping 


given  by 


X  =  y  = 


which  is  a  one-to-one  mapping  of  S  onto  ^ 

quadrilateral)  or  T  onto  ^  triangle.)  The  functions 

and  are  analytic  functions  in  tj  and  have  certain  uniform 

bounds  of  the  form 


(4.1a) 


|d“x^.'^^|.  |d“y^.^^1  s  C.h  ..  |a|  =  1.2 

'  J  '  '  J  '  1  n.j*  I  I 


(4.1b) 


C.-h^  .  <  IJ.I  <  C-h^  . 

2  n,j  '  j'  3  n.j 


where  Jj  is  the  Jacobian  and  the  constants  i  =  1,2,3,  are  independent 


40 


of  the  partition.  Here,  h  .  is  the  maximal  side  length  of  x. 

n.j  *  j 


(n) 


There 


are  some  additional  assumptions  about  uniform  analytic  extendability  of 

functions  and  outside  S  or  T,  which  we  will  not 

0  J 


(n) 


are  straight  sided 


discuss  here  (see  [BGl]).  In  the  case  when  r. 

<3 

triangles,  the  mapping  is  linear  and  the  uniformity  condition  is 

satisfied  if  the  minimal  (and  maximal)  angle  is  bounded  from  below  (and  above] 
by  constants  independent  of  the  partition  under  consideration.  Further,  we 


assume  that  x^.^^  n  x,^^^ 
J  k 


is  either  empty,  a  single  common  vertex  or  a  single 


-(n) 


common  side  and  if  s  .  ,  =  x . 

J.k  j 


-(n) 


is  the  common  side,  then  x  €  s .  , 

J.k 


is  mapped  by  the  mappings  (F^.^^)  ^  and  (F^^^)  ^  into  the  same  relative 
position  on  the  side  of  the  standard  square  or  triangle.  The  above 
assumptions  about  the  partitioning  are  analogous  of  the  ones  for  the 
partitioning  of  (0,1)  into  Ij  in  Section  2. 

Now  we  will  define  the  quasiuniform,  radical  and  geometric  meshes. 

1)  The  guasiuniform  mesh.  Here  {7^}  is  such  that 


h  . 

n,  J 


£  X, 


where  h  =  max  h  .  and  x  >  0  is  uniform  for  all  meshes  7  under 
n  n,  J 

consideration. 

2)  The  radical  mesh  with  power  3^1  (with  respect  to  ) .  If 
i  with  jf  =  1  -  ^ 


C.h  min,  ,  r  (x) 
J 


s  h  .  i  C  _h  max  ,  .  r  ^(  x ) 
n.J  Z-nt^^-In) 

J 


and  if  A.  € 

1  J 


then 


41 


'  h  I 

1  n[ 


max,  .  r^(x) 
-(n) 
xex  . 

J 


s  h  .  s  C_h  max,  .  r^Cx] 

n.J  Z-nt^^-Cn) 

j 


Here  r(x)  is  the  distance  of  x  from  This  definition,  which  holds  in 

a  neighborhood  of  A.,  is  the  two-dimensional  version  of  the  definition  of 

the  radical  mesh  introduced  in  Section  2.  It  was  employed  in  [BKP] . 

3)  The  geometric  mesh  with  ratio  q,  0  <  q  <  1  (in  the  neighborhood 

of  A.).  Let  us  number  the  elements  of  by  a  double  index  with  j  = 

1  j,k 

1 . p(k),  p(k)  s  p-  and  k  =  l,2,...,n+l.  Then  if  d  .  ,  denotes  the 

u  n ,  J  ,  K 

distance  between  and  A.  and  A.  i  ,  , 

J.k  1  1  n.j.k 

_  n+2-k  ^  ^  r-  n+l-k 

S'!  =  <*n,J.k  =  '=2" 

‘^l^^n.j.k  ■  ^n.j.k  “  '^z’^n.j.k 

J  =  1, . . . ,p(k),  k  =  2, . . . ,n+l. 

If  A.  €  then  k  =  1  and 

1  J.k 

^  \.J.k  ^  J  =  1..--.P(1). 

where  constants  C.  and  k.  are  uniform  for  all  meshes  under  consideration. 

1  1 

The  geometric  mesh  will  be  denoted  by  more  about  these  meshes  we 

refer  to  [BGl]. 

As  in  one  dimension,  the  radical  and  geometric  meshes  are  refined  in 
neighborhoods  of  the  vertices  A^^  where  the  singularities  of  the  solution  are 
located. 

Now  let  us  define 


(4.2)  V"(n)  =  {u  €  H^n)|  For  x  6  t[^\  u  o  €  U  (K)} 

J  J  Pj 
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where  K  s  T,  li  =  P  if  is  a  curvilinear  triangle  and  K  s  s, 

U  s  Q  or  Q'  if  and  is  a  curvilinear  quadrilateral.  Here, 

P  (K)  denotes  the  set  of  all  polynomials  of  total  degree  ^  p.  on  K  =  T 

or  K  =  S.  Further,  Q  (S)  is  the  set  of  all  polynomials  of  degree  s  p. 

P  •  P  • 

in  each  variable  and  Q'  (S)  =  P  (S)  ®  The  elements  Q  are 

sometimes  called  product  elements,  while  Q'  are  called  serendipity  or  trunc 

elements,  see  e.g.  [Ci],[SzB].  They  are  used  for  example,  in  the  programs 
MSC/PROBE,  FIESTA,  PEGASYS  and  STRIPE. 


The  theoretical  question  of  under  what  conditions  the  space  Qp^S)  or 
Qp(S)  is  preferable  has  not  yet  been  resolved  sufficiently.  For  the 
analysis  of  some  aspects  of  this  problem,  we  refer  to  [BE]. 

We  now  set  v”  =  n  H^.  In  general,  the  condition  u  e  in  the 

definition  of  restricts  the  functions  u(Fj'^^(x))  to  belong  to  a 

proper  subset  of  P  (T)  (or  Q  (S),  Q'  (S)).  This  is  necessary  to  enforce 

continuity  of  u  €  for  example,  when  p^  is  different  over  adjacent 

elements. 

As  in  Section  2,  we  now  deal  with  the  space  V^(3’^,2^)  c  ^^(12)  where 
is  the  degree  vector.  We  note  that  for  quasiuniform  meshes  and  uniform 
=  p^,  we  have 

-2  2 

N  s  h  ^p  . 
n  n  n 


which  also  holds  for  the  case  of  radical  meshes  and  uniform  p  .  For  the 
geometric  mesh  and  uniform  p^  we  have 


N 

n 


nPoP^ 
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Analogously  as  in  Section  2,  we  have  in  our  model  problem  (3.1),  B(u,u)  = 
|(u|(p  »  hence  with  e  -  u  -  u  as  before. 


n"E 


inf 


llu 


vll 


4.2.  Approximation  theorems 

Let  us  now  formulate  theorems  analogous  to  the  ones  in  Section  2  from 
which  the  error  bounds  for  the  finite  element  method  immediately  follow. 

Let  us  remark  that  when  comparing  the  performance  of  different  elements 
like  Qp  and  elements  (for  example),  the  rate  of  convergence  with 

respect  to  p  or  h  in  the  theorems  below  may  be  the  same,  but  the 
performance  may  be  different  (see  [BE]).  First,  for  the  h-p  version  with 
quasiuniform  meshes,  we  have  the  following  theorem. 


Theorem  4. 1  [BSl].  Let  us  consider  =  p  and  a  family  of  quasiuniform 
meshes  7^.  Then 

(4.3)  inf  ||u  -  v||  <  Ch^"^p"^^'^^||u|| 

^^yn  H^n)  H  (n) 

where  p  =  min(p+l,k)  and  C  is  a  constant  independent  of  h,p  and  u. 

From  the  above  theorem,  it  is  clear  that  for  the  p  version,  we  have  a 
- (k-1 ) 

rate  of  0(p  ).  Let  us  outline  how  this  rate  is  established — this  will 

allow  us  to  explain  some  key  techniques  used  in  p  version  approximation 
theory. 

Suppose  we  consider  a  mesh  of  triangles,  then  for  the  p  version, 
there  will  only  be  a  finite  number  of  elements,  this  number  being  independent 
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of  p,  l.e.  n.  Consider  the  triangular  element  x.  (with  h  .  <  1). 

The  first  step  is  to  extend  the  function  u  defined  on  open  unit 

square  R  satisfying  R  o  such  that  the  norm-preserving  property 


k  C||u|i 

H^(R) 


J 


holds. 


Next,  we  use  the  Tchebysheff  transformation,  which  transforms  algebraic 
polynomials  into  periodic  trigonometric  polynomials,  u  gets  transformed  into 
u  on  the  mapped  square  R  such  that  the  norm  is  once  again  preserved,  i.e. 


H'"(R) 


The  function  u,  of  course,  is  now  periodic  and 


therefore  can  be  approximated  by  Fourier  expansions.  For  each  p  =  1,2,.  . 
we  therefore  obtain  a  trigonometric  polynomial  w  satisfying 


1^  -  “pll  s  -  Cp"^‘''^^||u|| 

P  H^(R)  H*"(R) 


,  0  £  s  s  k. 


-^(*)  (n) 

Transforming  the  functions  u,  w  back,  we  get  a  e  T  (x.  j  for  each 

P  P  P  J 

p,  J,  such  that 


(4.4) 


11^  J )  II  ^ '  -(k-s)  ||-»|| 

ll""  “  p  i  s  (n)  -  k  (n) 

P  H  (x.  H  (x."M 

J  J 


Note  that  in  the  above,  we  used  the  fact  that  is  a  proper  subset  of  R, 

l.e.  the  distance  of  from  the  boundary  of  R  is  strictly  positive. 

The  above  construction  yields  a  separate  polynomial  for  every  triangle. 

The  next  step  is  to  adjust  these  polynomials  so  that  they  will  be  continuous 
on  n.  As  shown  in  [BS2],  for  the  case  k  >  3/2,  we  may  first  adjust 

so  that  it  coincides  at  the  vertices  with  the  function  u  and  estimate  (4.4) 

is  preserved.  Then  the  Jump  (discontinuity)  on  a  side  y  of  x^.^^  is  a 


polynomial  which  vanishes  at  the  ends  of  j  and  satisfies  (by  (4.4)  and  the 
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trace  theorem), 


(4.6) 


„l/2,  ,  -  '“U.  (n), 

"oo  ' 


The  final  step  is  to  extend  ^  to  a  polynomial  over  which 

coincides  with  on  and  vanishes  on  the  other  two  sides  of  ts  and  which 
has  a  suitably  bounded  norm.  In  the  original  proof  from  [BSK]  , 

this  was  accomplished  in  a  non-optimal  way,  giving  an  estimate 


(4.7) 


I  p 


hVt*"') 

J 


H  (y) 


The  right  side  of  (4.7)  was  bounded  above  by  Cp  ^^||u||  r  1  that 

H^(t.  ^) 

J 

subtracting  from  each  j  gave  a  polynomial  w^  € 

satisfying 


(4.8) 


II 

u  -  W 


-  Cp  u 


H^(n) 


H^(n) 


Hence,  we  obtained  the  desired  estimate,  except  that  a  factor  p  ^  was  lost. 

— c 

This  loss  can  be  reduced  to  p  ,  for  e  >  0  arbitrary,  by  a  useful 
interpolation  argument,  which  we  illustrate  below. 

We  obviously  have 


(4.9) 


^inf  ||u  -  ^\\  , 

weV"  H 


lu  -  P  u 


(n) 


D  "  1 

P  H^(n) 


s  Cllul 


H^(n) 


where  PpU  is  the  projection  of  u  onto  V^.  Also,  from  (4.8),  we  get 


(4. 10) 


-  PpUf  1 
P  HMfl) 


^  „  -(k-2),|  I 

i  Cp  u 


H‘'(n) 
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The  above  interpolation  argument  has  been  used  frequently  in  the  p  version 
literature,  e.g.  in  [BSK] ,  [Dl],  [JV],  [SI].  Note  that  such  an  argument  does 

not  work  for  the  h  version,  since  the  rate  of  convergence  is  bounded  by  the 
polynomial  degree. 

To  obtain  the  optimal  order  (4.3),  without  the  c,  the  extensions 
in  (4.7)  must  be  constructed  in  an  optimal  way.  This  was  accomplished  in 
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[BSl],  using  the  theorems  4.3  and  4.4  below,  which  have  found  applications  in 
other  contexts  as  well  (see  [BS3] ,  for  example).  Using  these  theorems 
together  with  (4.6)  gives  the  desired  optimal  estimate  for  leading  to 

(4.3). 

Theorem  4.3.  [BSl],  [BCM] .  Let  T  be  the  standard  triangle  with  sides 

i  =  l.?,3,  and  f  oe  coiitinuous  on  oT  with  f!  e  P  (y.),  i  =  1,2,3. 

pi 

Then  there  exists  U  e  P  (T)  such  that  U  =  f  on  3T  and 

P 


H  (T) 


I  1/2 
H^^^(aT) 


where  the  constant  C  is  independent  of  p  and  f. 


Theorem  4.4.  [BSl].  Let  S  be  the  standard  square  with  sides  i  = 

1,2, 3, 4.  Let  f  be  a  continuous  function  on  dS  with  f|  e  P  (y.), 

P  1 

i  =  1,2, 3, 4.  Then  there  exists  U  e  Q  (S)  such  that  U  =  f  on  3S  and 

P 


I  1 

H  (T) 


'  1/2 

H  (as)’ 


where  the  constant  C  is  independent  of  p  and  f.  (See  also  [BCM]) 
Theorems  4.3  and  4.4  are  the  analogs  of  the  classical  extension 
theorems. Here  the  extension  is  constrained  to  polynomials  of  the  same  degree 
as  the  traces.  We  mention  that  when  using  the  space  Q'p(S)  instead  of  0^(3), 
the  constant  C  in  Theorem  4.4  is  mesh  dependent  ([Md]). 

Let  us  now  consider  singular  functions. 

Theorem  4.5.  [BSl].  Let  u  =  r**  log^r  t/f(e)  where  by  r*^,  a  complex,  we 
denote  Re  r**  or  Im  r°^,  and  let  the  meshes  be  quasiuniform  with  uniform  p. 
Then 
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[  min(a,  p-a) 

h“,  - ^ -  . 

2a  J 

vev  IJ  ,pj  p 

where  g(h,p,s)  =  inax(|log  h|  ,  |  log  p]  ),  a  =  Re  a  and  C  is  a  constant 
independent  of  h.p. 

See  also  [BD] .  The  next  theorem  is  analogous  to  Theorem  2.5  and  gives 
the  rate  of  convergence  of  the  h  version  with  radical  meshes. 

Theorem  4.6.  Let  u  =  r^  log^r  i^(0).  Then  for  the  radical  mesh  with  power  B 
and  uniform  p, 

inf  ||u  -  v||  i  C(a,B,s)h^ 

,,n  H^(n) 

veV 

provided  that  l>y=i-i>i-^>o. 

The  proof  is  a  simple  generalization  of  the  results  in  [BKP] . 

Finally,  let  us  state  the  following  theorem  which  shows  that  the  h-p 
version  with  geometric  meshes  leads  to  exponential  convergence. 

Theorem  4.7  [BGl].  Let  u  €  ,  0  <  0  <  1.  Let  be  geometric 

P 

meshes  and  junspsi/n,  p2l,  0<p,i^<oo.  Then  we  have 

-3r?Tr 

(4.11)  inf  ||u  -  v|  s  Ce  ^ 

H  (Q) 

with  j  >  0. 

See  also  [BD] .  The  proof  of  Theorem  4.7  uses  a  combination  of  the 
approaches  used  in  the  outline  of  the  proof  of  the  convergence  of  the  p 
version  in  this  section  and  the  proof  in  the  one  dimensional  setting  presented 
in  Section  2  for  the  h-p  version  exponential  rate. 

It  is  more  advantageous  to  use  a  linear  degree  vector  than  a  constant 
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one.  The  use  of  the  geometric  mesh  is  advantageous  when  the  effort  of  mesh 
generation  is  considered  for  complex  geometries,  since  the  refinement  is 
restricted  to  a  small  subset.  It  is  also  possible  to  use  a  sequence  of 
radical  meshes  as  in  one  dimension  and  by  means  of  the  envelope  obtain  an 
exponential  rate,  but  the  mesh  generation  effort  is  now  much  more  complex. 

So  far  we  mentioned  only  the  problems  where  the  approximation  in  the 
spaces  H^(Q)  was  relevant.  There  are  problems  where  the  norms  in 

have  to  be  used.  Various  analogous  results  can  be  proven  here.  We 
refer  to  [Go]  [S2]  for  details.  An  analysis  of  the  error  in  the 
and  negative  norms  has  been  presented  in  [JS] .  Let  us  briefly  describe  a 
result  of  interest  from  this  reference,  related  to  the  duality  results 
presented  in  Section  2. 

Consider  the  L-shaped  domain  in  Figure  4.1.  The  solution  of  Laplace’s 

2/3 

equation  on  n  is  known  to  have  a  singularity  behaving  like  r  by  the 
results  of  the  previous  section.  If  the  h  version  is  used  to  approximate  the 


B  A 


Figure  4.1.  The  L-shaped  domain. 


solution  of  Laplace’s  equation  (using  a  uniform  mesh),  the  error  will  satisfy 
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The  usual  duality  argument  does  not  lead  to  a  full  extra  power  in  h  for  the 
L^(Q)  norm  over  the  H^(Q)  norm.  This  is  because  of  the  non-convexity  of 
the  domain;  the  solution  of  the  dual  problem 

(4. 13)  -Aw  =  e  on  Q 


will  satisfy 


(4. 14) 


k  -  ^  ^'il  (n) 


with  k  =  4/3  only  and  not  k  =  2,  as  required.  We  note  that  the  estimates 
in  (4.12)  cannot  be  improved,  see  (W1 ] . 

For  the  p  version,  however,  we  obtain  the  full  extra  power  of 


ronvergence,  i.e. 


(4. 15) 


1 

H^n) 


I^Il  (n)  “  1 

2  ^  H  (n) 


s  Cp 


This  is  because  the  solution  of  (4.13)  can  now  be  decomposed  into  a  singular 

and  smooth  part  as  in  equation  (3.10)'.  The  smooth  part  will  satisfy  (4.14) 

2/3 

with  k  =  2.  The  singular  part  will  not;  however,  since  it  behaves  like  r 
it  is  approximated  at  double  the  usual  rate  by  the  p  version,  and  the  duality 
argument  easily  follows.  The  same  argument  works  for  other  non-convex  domains 
( see  [ JS 1  '  . 

So  far,  we  have  assumed  that  g  =  0  on  r^.  In  the  inhomogeneous 
case,  u  is  approximated  by  first  finding  a  suitable  approximation  g^  to  g 
in  the  trace  space  w”  =  w'^(rj^)  =  {wj^.  ,  w  €  v"^}  where  is  the  finite 

element  space  when  =  <p.  This  introduces  a  new  component,  depending  upon 

the  difference  g  -  g^,  into  the  error.  There  are  various  interesting 
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results  on  the  effect  of  this  component  on  the  total  error,  depending  upon  how 
is  defined.  We  refer  to  [BS3]  for  example. 

Let  us  mention  here  that  the  spectral  element  method,  which  has  many 
features  in  common  with  the  p  version,  also  leads  to  various  approximation 
results.  We  will  not  survey  them  here  but  refer  instead  to  [MdP]  and  the 
references  therein. 

In  the  three-dimensional  case  the  solution  is  essentially  similar  but 
more  complex.  Let  us  mention  the  following  theorem  for  the  p  version. 


Theorem  4. 8  [D2].  Consider  a  fixed  mesh  J  and  a  sequence  of  degrees  p^  ^ 

0.  Let  f  €  C^Cn).  If  u  is  the  solution  of  (3.13)  and  u).,  .  ,  J  = 

J  (-,1 

l,...,n^,  I  =  l,...,n^,  i  =  1,2, .  are  as  described  in  Section  3,  then 


for  arbitrary  c  >  0 


where 


inf  ||u  -  v||  £  Cp 

veV  H  (fl) 


min{|^,  } 

j.«"j 


and  C  depends  on  u  and  c  but  is  independent  of  p. 

The  above  theorem  is  based  on  the  regularity  results  (3.14),  (3.15)  for 
the  solution  u. 

Theorem  4.5  indicates  what  was  observed  computationally  for  the  model 
problems  we  considered,  namely  that  in  most  practical  cases  the  edge 
singularity  (and  not  the  vertex  one)  governs  the  rate  of  convergence  for  the  p 
version. 

In  three  dimensions  the  proper  mesh  selection  is  complex  because  of  the 
(anisotropic)  singular  behavior  of  the  solution.  Here,  in  the  neighborhood  of 
edges,  long  (needle)  elements  have  to  be  used,  because  while  the  solution  is 
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smooth  along  the  edge,  it  is  unsmooth  in  the  direction  perpendicular  to  the 
edge.  In  the  neighborhood  of  the  vertices,  the  singularity  is  isotropic  and 
elements  of  the  usual  type  are  used.  Hence,  the  class  of  meshes  (radical  and 
geometric)  introduced  in  two  dimensions  has  to  be  enlarged.  In  addition, 
the  mesh  generation  has  to  be  more  sophisticated,  to  take  into  account  the 
transition  between  various  meshes.  For  the  approximation  properties  of 
"needle"  elements  which  violate  the  usual  minimal  angle  condition,  we  refer 
to  [BA],  [Krl. 

As  we  saw  in  Theorem  4.7,  for  a  proper  sequence  of  meshes  and  degree 
distributions,  we  have  an  exponential  rate  of  convergence  (with  the  h-p 
version).  In  the  three  dimensional  case,  we  get  a  result  which  is  quite 
analogous  to  Theorem  4.7,  with  the  exponent  v'FT  in  (4.11)  now  being  replaced 
by  .  For  more  see  [BG5]. 


4.3.  The  effect  of  numerical  quadrature 

So  far,  we  have  assumed  throughout  that  exact  integration  is  used  in  the 
calculations  of  stiffness  matrices  and  load  vectors.  In  practice,  of  course, 
these  are  generally  computed  using  numerical  integration.  Unlike  the  h 
version,  where  a  rule  of  fixed  precision  may  be  used  as  h  ^  0,  in  the  p 
version,  the  precision  of  the  rule  must  also  increase  as  p  oo.  In  [BnS]  , 
the  optimal  convergence  of  the  p  version  using  numerical  integration  in  , 
t  =  1,2,3  was  established  under  a  set  of  sufficient  conditions  on  the 
quadrature  rule.  Let  us  describe  these  results  for  the  model  problem  (3.6), 
where  s  F^,  g  =  0,  a  fixed  mesh  J  =  7  is  used  and  p^  =  p  in  (4.2). 

Suppose  K  is  a  reference  element,  K  =  S  or  T.  Then  a  quadrature 


rule  on  K  induces  a  corresponding  rule  on  each  element  We  then  use  a 

family  of  rules  satisfying  the  following  two  assumptions. 
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(A)  The  weights  are  positive  and  the  quadrature  points  lie  within  K. 

(B)  R  is  exact  for  all  v  e  11  (K)  with  m  z-  2p. 

p  m 

Condition  (B)  can  be  weakened,  see  [BnS] . 

Newton-Cotes  rules  are  unsuitable  for  the  p  version,  since  they  have  low 

degree  of  precision  and  violate  the  positivity  of  weights  when  p  is  large. 

The  most  advantageous  rules  are  Gaussian  rules,  which  obviously  satisfy  (A). 

In  [BnS],  it  is  shown  that  in  1  -  d,  the  p  point  Gauss-Legendre  and  p  -  1 

point  Gauss-Lobatto  rules  are  the  minimal  required  when  polynomials  of  degree 

p  are  used.  For  Q  type  elements  in  higher  dimensions,  the  corresponding 

P 

tensor  product  rules  using  p  +  1  points  in  each  direction  are  minimally 

required.  For  other  rules,  see  [BnS]. 

Let  A  =  diag(k.],  where  k.  =  k. (x)  are  the  coefficients  in  (3.6). 

1  11 

Also,  let  DF.^(x)  denote  the  Jacobian  matrix  of  the  inverse  of  F.  and 
J  J 

define  the  matrix  B"^  =  B'^(x)  by 


B'’  =  [b^J  =  J.{DF:")A(DF:") 

J  J  J 


-1  -1  T 

j^)A(DFj  )^' 

Then  the  following  theorem  is  proved 


Moreover,  let  B  =  max  |lb,'^„|| 

^■0  J.k.f 

in  [BnS] . 

Theorem  4.9.  Let  £2  c  r'"  .  Let  f  in  (3.3)  be  in  (Q) ,  s  >  t/2  and  let 

k  1 

u,  the  solution  of  the  model  problem  (3.6),  lie  in  H  (£2),  k  >  1.  Let  bj^^  € 
H  (K)  for  each  J,  k,  i,  with  d  >  t/2.  Let  e  =  e^  be  the  error  when  the 
finite  element  solution  u^  is  computed  using  rules  which  satisfy 

(A),  (B).  Then  for  any  0  s  q  <  m-p,  r  =  min(p, m-p-q ) , 


I  1  -  c{(m- 

H^Q)  V 


-(s--)  -(d-^) 

^  l|f||  s  ^d  2  Ill'll  1 

H  (£2)  *  H  (£2) 
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+  r 


-(k-1) 


-(d--) 

^d,2j 


l"i  k  } 


where  the  constant  C  is  independent  of  u,  m,  p  and  q. 

The  above  estimate  contains  three  terms.  The  first  term,  which  is 

-(s — ) 

2 

0((m-p)  )  arises  from  the  approximation  of  f,  the  second,  which  is 

-(d-i) 

0(q  )  arises  from  approximating  the  mappings  and  coefficients  k^, 

-(k-1 ) 

and  the  final  term,  0(r  ),  arises  from  the  solution  u. 

For  the  case  that  f,  F^  and  k^(x)  are  sufficiently  smooth,  the 

above  theorem  shows  that  if  m  ~  2p,  then  we  once  again  recover  the 
-(k-1 ) 

0(p  )  convergence  predicted  by  Theorem  4.1.  A  similar  result  holds  for 

the  case  that  the  solution  is  singular,  for  which  it  may  be  shown  that  the 
optimal  rate  of  0(p  may  be  recovered. 

If,  however,  d  is  small,  then  "overintegration"  (i.e.  taking  m  >  2p) 
may  be  needed  to  preserve  the  error  bound.  See  [BnS],  (Ki]  for  additional 
theorems  and  various  computational  results. 


4.4.  Numerical  examples 

The  theorems  above  are  asymptotic  in  nature.  To  verify  that  they  hold 
for  a  practical  range  of  parameters  as  well,  we  consider  a  computational 
example  involving  the  elasticity  problem  (3.7)  over  the  L-shaped  domain  shown 
in  Figure  4.1. 

We  specify  tractions  on  dQ  (with  =  ^,  t  =  0  and  sides  OA,  OE 
being  traction  free)  such  that  the  true  solution  u  =  (u^,U2)  is  given  by 

(4.16)  u,  =  r**[  (K-Q(a+1 )  )cos  a0  -  a  cos(a-2)0]  (Mode  1) 

1  zu 
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^  ((c+Q(a+n  )sin  a0  +  a  sin(a-2)0] 


(Mode  2) 


u 


2 


where  a  =  0.5444837  and  Q  =  0.5430756.  Here,  G  is  the  modulus  of 
rigidity  and  k  =  3  -  iv,  with  v  =  0.3.  Equation  (4.16)  shows  that  the 
solution  has  an  r°^  type  singularity  at  the  origin  0,  of  the  same  type  as 


in  the  decomposition  (3.10)  -  (3.11).  It  may  also  be  mentioned  that  the 
solution  lies  in  ^(fl)  for  arbitrary  e  >  0  and  in  the  Besov  space 

2, CO 


We  first  consider  the  FEM  approximation  of  this  problem  using  a  uniform 
mesh  with  square  elements,  as  shown  in  Figure  4.2. 

Since  we  have  a  quasiuniform  mesh,  by  Theorem  4.5,  for  p  ^  1  we  obtain 
the  estimate 


|e||^  s  C  min 


in  h  ,  — 


min(a,p-a) 


2a 


This  leads  to  the  predicted  asymptotic  rates  h°^  and  p  for  the  h  and 


p  versions,  respectively. 
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Figures  4.3  and  4.4  in  which  the  relative  error  in  the  energy  norm 
|e|g.  is  plotted  respectively  as  a  function  of  h  and  p,  show  that  the 
error  is  in  the  asymptotic  range  even  for  moderate  h  and  p.  We  reach  the 
same  conclusions  as  in  the  one-dimensional  case. 


2  4  6  8  10 


MESH  SIZE 

Figure  4.3.  The  error  for  the  h  version. 


DEGREE  p  OF  ELEMENTS 

Figure  4.4.  The  error  for  the  p  version. 
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DETAIL 


Figure  4.5.  The  geometric  mesh. 


Let  us  now  consider  the  case  of  strongly  refined  (geometric)  meshes. 

Such  a  mesh  is  shown  in  Figure  4.5  for  n  =  2  layers.  This  is  a  geometric 
mesh  with  ratio  0.15  (which  leads  to  nearly  optimal  convergence). 

Figure  4.6  shows  the  relative  error  vs.  the  number  of  degrees  of  freedom 
for  various  n  (with  p  being  increased  in  each  case  to  increase  N)  in 
log||e|j^  vs  N  scale.  In  Figure  4.7,  log  ||e||p.  ^  for  various  combina¬ 
tions  of  n  and  p  has  been  plotted.  It  is  observed  that  pj  ~ 

Ce  ,  which  is  consistent  with  Theorem  4.7. 

Remark  4.1.  The  meshing  required  by  the  h-p  version  can  be  expensive  to 
implement.  Often  (especially  in  the  context  of  commercial  p  version  codes), 
it  is  enough  to  select  a  fixed  "good"  mesh,  i.e.  one  that  is  sufficiently 
refined  near  potential  singularities,  so  that  the  erVor  decreases  largely  in 
the  pre-asymptotic  exponential  region.  This  region  is  clearly  seen  in  Fig. 

4.6  for  any  fixed  n  —  when  p  is  increased  further,  the  error  goes  into  the 
asymptotic  region,  which  is  algebraic.  In  general,  for  small  required 
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accuracies,  the  mesh  should  be  over-refined  (low  p),  while  for  high 
accuracies,  it  is  under-refined  (high  p).  The  optimal  design  of  the  mesh 
would  required  either  a  fully  adaptive  approach  (see  e.g.  [DO])  or  an  approach 
based  on  an  expert  system  (see  e.g.  [BR]). 

We  presented  here  only  examples  of  two  dimensional  problems.  For  some 
computational  results  of  three  dimensional  problems,  we  refer,  for  example,  to 
[AB] ,  [CM].  The  basic  features  are  completely  similar  to  the  two  dimensional 
case. 
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Figure  4.6.  The  error  )le|g.  ^  as  a  function  of  n  and  p. 
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Figure  4.7. 


The  error  ll^i^  selected  combinations  of 


(n,  p)  . 
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4.5.  The  pollution  problem 

In  Section  2,  we  discussed  briefly  the  "pollution  effect"  in  the  p 
version.  By  "pollution",  we  mean  the  phenomenon  by  which  a  singularity  in  the 
solution  at  one  point  influences  the  behavior  of  the  approximation  over  the 
entire  element.  With  large  elements,  the  accuracy  can  thereby  be  unfavorably 
affected  even  at  points  that  are  distant  from  the  singularity.  This  effect  is 
even  more  pronounced  when  it  is  the  derivatives  of  the  solution  (like 
stresses)  that  are  of  interest.  Some  mathematical  analysis  of  pollution  may 
be  found  in  [Wl],  however,  a  complete  theory  is  not  yet  available. 

Let  us  show  the  effect  of  this  pollution  by  an  example  from  [BO]. 
Consider  the  plane  strain  problem  (E  =  10^,  v  =  0.3)  of  a  plane  loaded  by  a 
concentrated  load  as  shown  in  Figure  4.8.  The  boundary  conditions  are 
traction  conditions  such  that  the  solution  is  the  well  known  Bousinesgue 
solution  on  a  half  plane.  The  solution  then  has  infinite  energy  but  it  has 
finite  energy  on  any  subdomain,  particularly  in  the  shaded  area  shown  in 
Figure  4.8. 

In  the  case  of  the  h  version,  the  energy  error  in  the  shaded  area  will 
be  essentially  the  same  as  the  error  of  the  best  approximation  on  the  domain 
(let  us  call  it  pollution  free).  This  property  of  the  h  version  is  caused  by 
the  filtering  out  of  the  singularity  through  a  layer  of  elements.  The  same 
effect  of  filtering  occurs  in  the  h-p  version,  as  long  as  there  are  a 
sufficient  number  of  layers  of  elements  between  the  singularity  and  the  shaded 
domain.  In  Figure  4.9,  we  show  a  sequence  of  four  meshes  (on  half  of  the 
panel  because  of  symmetry)  used  for  the  computation  of  the  error  in  the  shaded 
area.  Figure  4.10  shows  the  relative  error.  We  have  also  depicted  the 
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pollution-free  error  which  is  the  error  of  the  best  approximation  of  the 


solution  on  the  shaded  area.  We  see  clearly  the  effect  of  the  pollution, 
particularly  when  the  p  version  is  used  on  an  insufficiently  refined  mesh,  as 
well  as  the  effect  on  this  pollution  when  the  layers  are  increased. 


P-DECxREE 

Figure  4.10.  The  error  in  the  shaded  area  for 
different  numbers  of  layers. 


5.  LOCKING  AND  ROBUSTNESS 

Suppose  we  consider  (3.7)  when  the  Poisson  ratio  v  is  very  close  to 
0.5  and  use  the  h  version  with  piecewise  linear  elements  to  approximate  this 
problem.  Then  it  is  well-known  that  the  finite  element  solution  for  practical 
choices  of  the  parameter  h  is  generally  not  very  accurate.  This  is  due  to  a 
phenomenon  called  locking  (Poisson  ratio  locking  in  the  particular  case 
mentioned),  which  is  said  to  occur  when  the  accuracy  of  numerical  schemes  for 
the  approximation  of  certain  parameter-dependent  problems  deteriorates  when 
the  parameter  is  close  to  a  limiting  value.  Another  example  of  locking  occurs 
when  the  h  version  is  used  on  (3.6)  when  the  parameter  d  =  is  closd  to 

0  (highly  anisotropic  materials).  Other  examples  of  locking  include  various 
plate  and  shell  models  when  the  thickness  t  is  close  to  0  (see,  for 
example,  [AF2],  [Pt],  [BrF]). 
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For  such  problems,  one  needs  methods  that  are  robust .  i.e.  which  work 


equally  well,  essentially  Independent  of  the  parameter  value.  In  this 
connection,  it  has  been  demonstrated  for  various  problems  that  the  p  version 
is  often  robust  in  situations  where  low-order  standard  h  version  schemes  are 
not.  For  example,  in  [Vo],  it  has  been  shown  that  in  the  case  of  Poisson 
ratio  locking,  the  p  version  is  essentially  free  from  locking  for  the 
displacements,  while  in  [Li],  the  Timoshenko  beam  model  (for  which  the 
standard  h  version  shows  locking,  see  [Ar] )  was  shown  to  be  free  from 
locking  when  the  p  and  h-p  versions  are  used.  Analogous  results  have  been 
shown  for  shell  models,  see  [Pt]. 

Let  us  point  out  that  other  robust,  locking-free  methods  exist.  For 
example,  in  [SV],  it  was  established  that  for  the  case  of  Poisson  ratio 
locking,  using  the  h  version  with  high  order  elements  (p  2  4)  also  avoids 
locking.  Another  strategy  to  overcome  locking  is  to  use  mixed  methods,  as 
done  for  the  Timoshenko  beam  problem  in  [Ar]  and  for  various  other  problems  in 
e.g.  [AF2]  [Pt]  [BrF].  However,  the  reformulations  of  the  variational  form 
necessary  for  mixed  methods  may  not  be  feasible  in  the  context  of  available 
commercial  codes,  so  that  it  may  be  more  desirable  to  have  a  method  (like  the 
p  version)  which  is  robust  in  terms  of  th^  standard  formulation  of  the 
problem. 

In  [BS4] ,  we  have  presented  a  general  theory  of  locking  and  robustness. 
Let  us  present  some  definitions  and  results  from  there,  adapted  to  the  context 
of  the  model  problem  (3.7),  when  the  material  is  nearly  incompressible. 

Let  us  denote  the  bilinear  form  (3.7)  by  B  (♦,•)  and  let  the  exact 

V 

solution  be  u  when  the  Poisson  ratio  is  v  (0  £  v  <  0.5).  Then  the 
corresponding  finite  element  solution  u|^  will  satisfy 
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B  =  B  (a'^.v) 

1/  n  If 


Vv  e  V 


Let  us  denote  the  energy  norm  (B^(w,w))^^^  by  ||w||g.  We  will  restrict  our 

attention  to  the  case  when  the  exact  solutions  u*^  belong  to  the  closed  unit 
B  k 

balls  H,  c  H  (n)  (k  i  2)  defined  in  terms  of  the  norms 
k,  v 


:|w||^  ,  =  l|wi\ 


+  {l-2v)  ^||div  wj|^ 

H  (n) 


As  shown  in  [BS5] ,  ^  correspond  to  precisely  the  correct  weighted  spaces 

that  characterize  the  regularity  of  the  solution  when  the  data  is 
appropriately  bounded. 

We  assume  that  the  sequence  {V^}  is  .such  that 


(5.2) 


C^F„(N^)  i  sup  inf  Jw  -  v||  5 

wsi}“  Ssv" 
k 


B  k 

where  H“  is  the  unit  ball  in  the  H  (Q)  norm.  Here,  are 

independent  of  N  and  Fq(N)  0  as  N  oo.  Equations  (5.1),  (5.2)  imply 

that  for  any  0  <  u'q  <  0.5,  we  can  find  constants  that 

the  following  holds  uniformly  for  all  0  s  i/  <  (N  =  N^), 


(5.3) 


C,(u„)Fo(N)  =  sup  s  C^IUolF^iN), 

k.i/ 


If  (5.3)  can  be  made  to  hold  uniformly  for  all  0^  <  0.5,  we  say  that  {V  > 

is  f ree  f rom  locking.  More  precisely,  we  have  the  following  definitions. 


Definition  5.1.  Let  f  be  an  increasing  function  with  lim  f(N)  =  oo.  The 

N-xx) 

extension  process  {V^>  shows  locking  of  order  f(N)  for  the  family  of 


problems  (5.1),  v  e  [0,0.5),  with  respect  to  the  solution  sets  in 
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the  energy  norm  iff 


<  lim  sup 

/ 

sup 

sup 

■1“  -  “JeJ 

N-«o 

i'6[0,0.5) 

> 

(FJNiflN)  ) 


-1 


=  M  <  00. 


k.i' 


For  the  case  that  M  (depending  on  f )  is  bounded  (respectively,  infinite), 
we  say  that  the  order  of  locking  is  ad  most  (respectively,  least )  f(N). 

If  (5.4)  holds  with  f(N)  h  l,  we  say  that  {V^}  is  f ree  from  locking. 


Def ini t ion  5.2.  Let  g  be  a  decreasing  function,  g(N)  0  as  N  ^  oo.  The 
extension  process  {V^}  is  robust  with  uniform  order  g(N)  for  the  family  of 

g 

problems  (5.1),  e  [0,0.5),  with  respect  to  the  solution  sets  in 


the  energy  norm  iff 


:5. 5) 


1  im 
N-<o 


sup 


sup 

U€[0,0. 


5) 


sup 

u  €H, 
k,  V 


I 

u 


-  U  il 

n'' 


(g(N) 


-1 


=  K  < 


We  see  from  the  above  definitions  that  if  f(N)  is  such  that  f{N)FQ(N) 
=  g(N)  -^0  as  N  00,  then  {V*^>  shows  locking  or  order  f(N)  iff  it  is 
robust  with  maximum  uniform  order  g(N).  Note  that  in  the  above  definitions, 
we  could  choose  other  error  measures  instead  of  the  energy  norm  (e.g.  the 
H^(Q)  norm),  see  [BS41. 

As  u-  0.5,  we  see  that  the  incompress ibi  1  i  ty  constraint . 

(5.6)  Cw=divw=0 


gets  imposed  on  the  exact  solution.  The  reason  locking  occurs  is  that  (5.6) 
gets  imposed  on  the  finite  element  solution  as  well.  For  there  to  be  no 
locking  in  the  energy  or  H^(Q)  norm,  the  following  necessary  condition  must 
be  satisfied 
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(5.7) 


CiFqCN] 


s  sup 


inf 


II  u  -  v|| 


->  ,,B  ->  ,,n 

ueH,  ^  veV 
k,  0 

Cv=0 


^  C^FqCN) 


where  „  =  n  H?  .  Condition  (5.7)  is  sufficient  as  well,  under 

certain  general  conditions  (see  [BS4] )  which  are  shown  in  [BS5]  to  be 
satisfied  for  our  model  problem.  This  reduces  the  whole  question  of  locking 
to  one  of  approximabi li ty  alone.  Let  us  give  some  sample  results  from  [BS5] 
in  this  connection. 

Let  us  take  Q  to  be  a  square  domain.  Let  and  be  the  uniform 

triangular  and  uniform  square  meshes,  respectively  (Figure  5.1).  Define,  for 


i  =  1.2,3, 


/^  ,  =  {u  €  C°(n),  ul  €  R^(S),  V  S  c 

p,  h  '  S  p  1 


where  R^(S)  =  (P  (S),  R^(S)  =  Q  (S),  R^(S)  =  Q' (S)  and  7!?  =  j!?. 

PPPPPP  32 


Figure  5.1.  Uniform  triangular  and  square  mesh. 
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Theorem  5.1.  Let  the  extension  process  be  the  h  version  for  the  problem 
(3.7),  (5.1)  using  subspaces  =  (V^  ^)^  with  p  fixed.  Let  the  solution 

g 

sets  be  ^  and  let  k  i  p  +  1.  Then  the  following  is  true  for  locking  in 
the  energy  norm  and  also  in  the  norm. 


Space 

Degree  p 

Order  of  locking  i 
f(N)  ~  CN^ 

Robustness  ; 

order  r 
g(N)  ~  CN“^  I 

(V^  )^ 

p.h 

1  <  p  <  3 

p  i  4 

1  =  1/2 

£  =  0 

r  =  (p-l)/2  ! 

r  =  p/2 

2  2 
(V^  .  ) 
p,h 

P  i  1 

£  =  1/2 

r  =  (p-l)/2 

3  2 

(V  ^) 

P.h 

1  5  p  <  2 

p  2:  3 

£  =  1/2 

£  =  1 

r  =  (p-1 )/2 

r  =  (p-2}/2 

[ 

The  above  shows  that  locking  cannot  be  avoided  for  p  £  3  for  the  h 

version.  Also,  using  rectangular  elements  leads  to  locking  for  all  p.  The 

1  2 

result  for  p  ^  4  for  (V^  ^)  was  first  proved  in  [SV]  and  holds  for  more 
general  meshes. 

In  contiast,  as  mentioned  earlier,  it  was  shown  in  [Vo]  that  the  p 
version  (using  straight  sided  triangles)  leads  to  uniform  robustness  of 
optimal  order  for  the  displacements,  so  that  locking  is  eliminated.  However, 
the  use  of  curvilinear  elements  is  indispensable  for  the  p  version  in  most 
cases.  In  this  connection,  we  have  the  following  computational  results. 

Let  us  consider  the  p  version  for  a  single  curvilinear  element  of  type 
Qp  for  problem  (3.7).  We  assume  that  this  element  is  the  image  of  the 
standard  square  under  a  curvilinear  mapping  F.  Below,  we  show  the  relative 


errors  in  the  energy  norm  for  the  cases  that  the  mapped  element  is  a  square,  a 


trapezoid,  a  quadrilateral  with  one  side  a  parabola,  and  a  quadrilateral  with 

one  side  a  circle.  We  do  this  for  u  =  0.3  (Figure  5.2)  and  i>  =  0.5  -  10 

(Figure  5.3).  It  is  observed  that  the  rate  of  convergence  now  behaves  like 
- (k-1 ) 

(p-s)  ,  with  s  s  0  depending  on  F  (s  =  0  for  affine  F).  There  is 

no  change  in  the  apparent  asymptotic  rate  so  that  in  this  sense,  there  is  no 
observed  locking. 

We  have  treated  various  other  examples  of  locking  in  [BS4]  using  the 
above  definitions  and  principles.  Like  the  above  case,  the  reason  locking 
occurs  in  general  is  due  to  the  imposition  of  a  constraint  similar  to  (5.6) 
Such  constraints  usually  can  be  interpreted  as  inter-element  continuity  re¬ 
quirements,  with  (5.7)  breaking  down  when  the  constrained  finite  element  space 
is  too  small.  (For  example,  (5.6)  can  be  interpreted  as  imposing 
continuity  on  a  space  related  to  V^,  see  [BS5].)  Unlike  low-order  h  version 
schemes,  the  p  version  remains  robust  under  such  continuity  constraints,  which 
is  why  it  is  more  resistant  to  locking. 


As  another  example,  consider  problem  (3.6)  with  ^^^^2  zero. 


In  this  case,  the  constraint  is  -5 —  =  0.  Suppose  the  h  version  is  used,  with 
piecewise  polynomials  of  degree  s  p  on  either  of  the  meshes  in  Figure  5.1. 


NUMBER  OF  DEGREES  OF  FREEDOM 
Fig.  5.2.  Error  behavior  for  curved  elements  for  u  =  0.3. 
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Fig.  5.4.  Rotated  triangular  and  square  meshes. 


Then,  since  the  meshes  are  aligned  with  the  directions,  there  is  no 

locking  (see  [BS4] )  since  there  are  a  number  of  functions  satisfying  ^  -  0. 

However,  if  these  meshes  are  not  aligned  with  x^  -  x^  direction  (for 
example,  as  in  Figure  5.4),  then  the  only  functions  satisfying  =  0 

be  actual  (rather  than  piecewise)  polynomials  and  (5.7)  will  not  hold.  There 
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will  be  complete  locking  in  this  case,  with  the  method  being  non-robust.  On 


the  other  hand,  if  we  use  the  p  version  with  just  one  element,  then  the 
dv 

constraint  ^ —  =  0  has  no  effect  and  there  is  no  locking,  i.e.  the  method  is 

robust.  See  [BS4]  for  a  detailed  analysis.  Also,  see  [SBS]  for  an  analysis 
of  locking  in  the  Reissner-Mindl in  plate. 


6.  EXTRACTION  AND  POST-PROCESSING  TECHNIQUES 

Usually  in  computational  practice,  the  solution  u  of  the  differential 
equation  considered  is  only  a  tool  to  get  the  primary  quantity  of  interest. 

For  example,  the  goal  of  the  computation  may  be  to  find  the  stresses  at  a 
point,  the  stress  intensity  factors,  or  the  resultant  (e.g.  the  moments)  in 
shell  and  plate  theory. 

The  results  in  Sections  2,4  indicate  the  high  asymptotic  rates  of 
convergence,  including  exponential  convergence,  that  can  be  obtained  by  the 
h-p  version,  when  the  error  is  measured  in  the  energy  norm.  The  question  that 
we  ask  here  is  how  these  rates  can  be  obtained  for  quantities  of  interest 
other  than  the  energy.  For  instance,  if  high-order  derivatives  of  the  h-p 
solution  were  computed  directly,  then  severe  oscillations  could  result  and 
these  rates  would  be  lost. 

Mathematically,  we  are  interested  in  evaluating  the  value  of  a  certain 
function  of  the  solution.  There  are  various  forms  in  which  this  function  may 
be  expressed,  using  properties  of  the  equation  under  consideration.  Choosing 
one  such  form,  we  then  substitute  the  finite  element  solution  instead  of  the 
exact  one  to  get  an  approximation  of  the  desired  value.  If  the  exact  solution 
were  substituted,  then  the  computed  value  would  be  the  same,  no  matter  which 
form  of  the  function  were  used.  However,  when  the  finite  element  method  is 
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used,  the  accuracy  of  the  computed  value  can  be  very  different.  In  this 
section,  we  indicate  how  the  convergence  rate  in  the  energy  can  be  achieved 
for  other  quantities  as  well. 

We  illustrate  the  main  idea  by  considering  the  simple  one-dimensional 
problem, 

-u"  +  u  =  f  on  Q  =  (0, 1 )  =  I 

(6.1) 

u(0)  =  u(l)  =  0 

2 

where  f  is  assumed  to  be  in  so  that  u  €  H  (Q).  Suppose  the 

quantity  we  are  interested  in  is  the  value  of  u'(0). 

Problem  (6.1)  can  be  cast  in  the  variational  form 


.1 


(6.21 


B(u,v)  = 


(u*v'  +  uv)dx  = 


fvdx 


with  u  €  Hj(n)  =  hJ.JI),  V  6  Hj|(n). 

Let  0  €  H^(£2j  satisfy  tpiO)  =  1,  0(1)  =  0.  Then 


B(u,0)  =  r  (-u"+  u)0dx  -  u'  (0) 

''O 


=  i! 


f0dx  -  u' (0) 


so  that 


(6.3) 


u'  (0)  =  r  f0dx  -  B(u,0) .. 

Jn 


Hence,  we  have  two  possible  ways  of  obtaining  the  desired  value  u'(0); 
we  could  either  use  differentiation,  or  equation  (6.3).  Function  0  is 
called  an  extraction  f unction.  Let  us  compare  the  accuracy  of  these  two 
approaches. 

Suppose  that  {V^>  is  the  sequence  of  subspaces  defining  the  extension 
process  (V*^  c  H^)  and  the  corresponding  finite  element 
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solutions.  Define 


(6.4) 


F  =  F(i 


u)  =  f  fij/dx  -  B(u,i/»)  (=  u'(0)) 
Jn 


(6.5) 


F  =  F( 
n 


u  )  =  f  f0dx  -  B(u  ,0) . 

n  Jq 


Then 


(6.6) 


IF  -  F  I  =  B(u  -u,i/»). 
'  n'  n 


Now  let  z  e  Hq(£2)  be  such  that 


(6.7] 


B(v,z)  =  B(v,^),  V  V  6  Hp(n) 


Then  for  any  e  we  have 


IF  -  F 


Hence , 


|B(u  -u,z)|  =  B(u  -u,z-w^^^)  s  C|lu  -  u||  ijz  -  , 

in'  n.  ''  n  "  1  '  'I  1 

H  («)  H^(n) 


|F  -  F  I  s  C||u  -  u||  inf  |z  -  a>|| 

H^Q)  wev"  H^n) 


and  the  theory  of  Section  2  is  applicable. 

The  smoothness  of  the  function  z  can  be  easily  established.  Usually, 

inf  ||z  -  u||  ~  ||u  -  u||  We  see  therefore  that  the  error  |F  -  F  | 

^  H  (n)  "  H  (n)  ^ 

weV^ 

2 

is,  in  general,  of  the  order  ||u  -  u|  ,  ,  i.e.  the  order  of  the  error  in 

^  hVq) 

the  energy.  In  contrast,  if  we  defined  F  by  direct  differentiation,  i.e. 

n 

F  =  u'(0),  then  we  would  only  get  ||u  -  u||  (modulo  a  logarithmic 

n 

term)  as  the  rate  of  convergence  of  jF  -  F  |. 

We  therefore  see  that  by  methods  of  the  type  above,  called  extraction 
methods .  more  accurate  values  can  be  obtained  when  0  is  properly  selected 
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ill)  smooth).  (In  the  case  above,  we  could,  for  instance,  take  0  =  1  -  x). 

Let  us  note  that  the  extracted  value  is  not  changed  when  we  replace  0  by  0 
-01,  we  V^,  so  that,  for  example,  in  the  h  version,  0  can  be  restricted 
only  to  the  element  containing  x  =  0.  (See  [BM]  for  more) 

Similar  ideas  can  be  used  for  other  values  of  interest.  One  important 
case  is  the  extraction  of  stress  intensity  factors  (see  Section  3).  Let  us 
elaborate  briefly  on  this  case.  Assume  that  we  are  solving  the  problem 

-Au  =  0  on  n 

where  Q  is  the  L-shaped  domain  shown  in  Figure  4.1.  Suppose  on  5Q  we  have 
the  boundary  condition 

=  g  on  ABODE 

on 

u  =  0  on  0E,0A. 

Then  using  (3.10)'  we  have 

_  2/3  ,20, 

u  =  F  r  cos(^)  +  V 

with 

I  I  I  -2/3,  I  .  I  I  1/3.  ,,  , 

I  V I  I r  I ,  I  grad  v| | r  |  =  o( 1 ) . 

Here,  F  is  the  stress  intensity  factor.  Let  us  define 

,  ,  -2/3  ,20,,  ,  , 

0  =  (r  cos(^)  ):^;(r ) 

where  a;  =  1  for  0sr<a,  ;^:  =  0  for  r>b,  a<b<l,  ;);e  c“(Q). 
Obviously,  w  =  A0  €  C^dl)  and 

w  =  0  for  r<a  and  r  >  b. 

Now  denoting  n  =  n  -  R  where  R  =  {(r,0)|  r  ^  c}  we  have 

C  C  E  ' 
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xiAifidx 


Jn  3n 


Hence,  taking  the  limit  as  c  ->  0, 


=  -(f  uwd.]->. 


Proceeding  as  before,  we  define 


F  =-ffuwdx]- 

n  [  in  ^  ^ 


and  get 


|F  -  F  I  d  C||u  -  u  II  Id  -  d  I 

H  (0)  H  (D) 


where  z  is  the  solution  of  the  problem 


-Az  =  w 


^  =  0  on  ABODE 

dn 


z  =  0  on  0E,0A. 


Because  w  €  and  all  internal  angles  except  at  0  are  z  is 

c"  except  in  a  neighborhood  of  the  origin.  In  this  case  (when  F  *  0)  we 


If  -  z  II  ^  ^  c|(u  -  u  |j 

^  H  (n)  H  (n) 

and  hence 


F  -  F 


n“  1 
H  (Q) 


(or  better),  i.e.  |F  -  F  |  is  once  more  of  the  order  of  the  error  in  the 

'  n ' 

energy.  This  is  typical  for  this  kind  of  extraction. 


The  stress  intensity  factor  extraction  method  obviously  requires  a 
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priori  knowledge  of  the  form  of  the  singularities.  In  two  dimensions,  a 
numerical  algorithm  for  obtaining  the  singular  form  for  the  general  elasticity 
problem  may  be  found  in  [Pa]  (see  also  [Sa]). 

As  a  numerical  example,  we  consider  the  L-shaped  domain  shown  in  Fig. 

4.1.  We  consider  the  elasticity  problem  described  in  Section  4  with 
prescribed  traction  conditions  on  Sn  such  that  the  exact  solution  is  given 
by 

u  =  Ui  ^  2u2 

where  u^,  u^  are  given  by  (4.16).  The  mesh  shown  in  Fig.  4.5  with  two 
layers  is  used.  In  Table  6.1,  we  report  results  using  the  p-extension  for  the 
given  mesh,  showing  the  number  of  degrees  of  freedom  N,  the  normalized  strain 


TABLE  1.  Strain  energy  and  stress  intensity  factors. 


p 

N 

STRAIN  ENERGY 

STRESS  INTENSITY 
FACTOR-MODE  1 

STRESS  INTENSTY 
FACTOR-MODE  2 

1 

41 

6.42072796 

0.95268 

2.29075 

2 

119 

6.74137580 

1.02177 

2.08422 

3 

209 

6.77029847 

1 . 00250 

2.02239 

4 

335 

6.77575144 

1 . 00073 

2.00437 

5 

497 

6.77683967 

0.99991 

2.00097 

6 

695 

6.77719530 

0.99985 

2.00022 

7 

929 

6.77736281 

0.99987 

2.00005 

8 

1199 

6.77749228 

0.99990 

2.00001 

00 

00 

6.77776914 

_ 

1 . 00000 

2.00000 

well  as  the  (normalized)  values  of  the  first  two  stress  intensity  factors.  In 
addition,  we  present  the  exact  values  (p  =  oo)  of  these  stress  intensity 
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factors.  In  Fig.  6.1,  we  show  the  relative  error  in  the  stress  intensity 
factors  as  well  as  the  relative  error  in  the  strain  energy  of  the  solution 
(i.e.  square  of  the  error  in  the  energy  norm).  We  see  that  the  errors  in  the 
stress  intensity  factors  are  in  fact  of  the  order  of  those  in  the  strain 
energy.  Note  the  typical  form  in  the  first  phase,  where  the  rate  is 
exponential,  and  in  the  second  phase,  where  it  is  algebraic.  For  more,  we 
refer  to  [SzB] .  For  the  extraction  of  stress  intensity  factors  in  3-d,  we 
refer  to  [AB] ,  [BPA] . 


'Z 


o 

Pi 

Pi 

W 

w 
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Fig.  6.1.  Error  of  the  stress  intensity  factors  and  strain  energy. 
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7. 


PLATE  MODELLING 


A  very  common  problem  in  structural  mechanics  is  the  problem  of  solving 
partial  differential  equations  on  a  thin  domain.  A  typical  problem  of  this 
type  is  the  plate  problem,  which  we  elaborate  on  below. 

Let  n  =  { (x  =  ,x 
w  is  a  bounded  domain  with  piecewise  smooth  boundary.  Denote 

^  WO  w 

S  =  {x  €  R  I  (x^  .x^)  e  d(j),  jx^l  <  R^=  {x  e  R  |  (x^  ,x^  )  e  w,  x^  =  ±  ^  K 

We  assume  that  d  <<  diam(u). 

We  are  interested  in  the  three-dimensional  elasticity  problem  which  we 
mentioned  in  Section  3.  Denoting 


2’^3^ 


(X1.X2) 


where 


and 

(7.1) 


u  =  (u^.u^.u^),  H^(n)  =  {u  e  (H^(n))^  |  u^  =  0  on  S} 


F(^)  =  g(x^,X2)(v3(x^.X2.5)  ^  V3(x^,X2.-|)dx^dX2, 


the  solution  u  €  of  the  plate  problem  satisfies  (3.4)  for  every 

V  e  HQ(n),  with  B(u,v)  given  by  (3.7)  (with  t  =  3).  This  problem  is  the 
so-called  plate  problem  with  soft  simple  support.  Using  other  constrained 
spaces,  we  can  describe  other  physical  conditions  (for  example,  using  the 
condition  u  =  0  instead  of  U3  =  0  on  S  gives  the  clamped  boundary 
condition) . 

In  order  to  numerically  approximate  this  problem,  a  common  technique  is 
to  first  replace  the  three-dimensional  model  by  a  two-dimensional  one,  which 
leads  to  a  considerable  reduction  in  the  number  of  degrees  of  freedom  and 
makes  the  problem  more  numerically  tractable.  The  simplest  model  that  could 
be  considered  is  that  of  the  Kirchoff  plate.  However,  this  has  the 
disadvantage  of  requiring  elements. 
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Two-dimensional  models  that  only  require  C  continuous  elements  can 
often  be  viewed  as  one  of  a  hierarchy  of  models  in  the  sense  described  below. 
To  obtain  such  models,  the  solution  is  essentially  "expanded"  in  terms  of 
polynomials  (or  other  functions)  in  the  x^  direction.  More  precisely,  let 
n  =  (n^.n^.n^),  n^^  s  0  integer  and  define 

^i  y2X 

(7.2)  "Hp(n)  =  {  u  e  HQ(n)  !  u.  =  ^ j  ^ 

j=0 

“3  ■  I  <'j'x,.X2>  ‘-j(V)} 

J=0 

where  Legendre  polynomials.  Then  by  the  hierarchy  of  plate  models 

we  will  understand  the  problem  of  finding  u^^^  such  that 

(7.3)  B(u^'^\v)  =  F(v)  V  V  €  ''Hj(fi). 

Let  us  note  that  because  of  symmetry,  the  models  (Zn^^  +  l ,  Zn^+l ,  2n^ )  are  the 
same  as  the  models  (2n^ +2 , 2n2+2, Zn^+l ) .  Typical  models  used  are  (n,n,n+l) 
models.  As  we  note  below,  such  hierarchies  of  models  are  easily  implemented 
in  the  framework  of  the  p  version.  The  exact  plate  solution  and  the  solution 
of  the  particular  model  depends  on  the  thickness  d.  It  is  possible  to  show 
that  the  models  n  =  (n^ , n^, n^) ,  n^  i  1 ,  n^  -  1  and  n^  i  2  converge  to  a 
limiting  solution  as  d  ->  0  which  is  the  same  s  the  limiting  solution  as  d 
0  of  the  exact  (three  dimensional)  solution  introduced  above.  The  model 
(1,1,0)  which  is  very  often  used  in  engineering  (the  so-called  Reissner- 
Mindlin  model)  does  not  have  the  above  property,  namely,  the  same  limiting 
solutions  as  the  exact  one,  unless  modified  elastic  constants  are  used  (see, 
e.g.  [BLl]).  We  note  that  for  v  =  0,  no  modification  is  needed.  Further,  it 
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used  in  (7.2)  is 


is  possible  to  show  that  the  polynomial  choice  in 
asymptotically  optimal. 

For  fixed  d,  the  solutions  of  the  models  converge  to  the  exact  3-d 
solution  as  min(n^ , n^, n^)  If  g(x)  in  (7.1)  satisfies  some  smoothness 

conditions,  then,  for  larger  n,  n  fixed,  the  rate  with  respect  to  d  (as 
d  ^  0)  is  higher.  The  various  models  have  different  properties  with  respect 
to  boundary  layers,  corner  singularities,  etc. ,  see  [BL2] . 

Assume  now  that  the  domain  u>  is  partitioned  into  the  elements  x. 

J 

based  on  triangles  T  or  squares  S  (see  Section  4).  Then  T^  =  x(-d/2,d/2)  is 

a  three-dimensional  element.  Now  denote  Qj(p,q)  to  be  the  set  of 
polynomials  of  degree  q  in  x  and  of  degree  p  in  x,,x  (this  could  include 

s3  i  c. 

the  case  of  P  Q  or  Q'  type  elements  in  x,,x_  —  see  Section  4).  Then 
the  finite  element  solution  with  elements  Qj(p,q)  of  the  three-dimensional 
plate  problem  with  degree  q  fixed  can  be  interpreted  as  the  finite  element 
solution  (using  elements  of  degree  p)  of  the  plate  model  with  n^  =  n.^  =  n^=  q 
in  (7.2),  Hence,  the  hierarchic  modelling  can  be  easily  implemented  in  the 
frame  of  the  p  (or  h-p)  version  of  tnree-dimensional  elasticity.  (In  fact, 
this  implementation  is  available  in  the  commercial  program  MSC/PROBE).  Let 
us  comment  on  three  main  aspects  of  the  solution  of  the  plate  problem  in  this 
hierarchical  framework. 

a)  The  locking  problem.  Here,  the  theory  addressed  in  Section  5  fully 
applies.  The  analog  to  constraint  (5.6)  for  the  (1,1,0)  model  is 


i,  1 


3^0 

ax. ' 

1 


i  =  1,2. 


For  more,  see  [SBS] . 

b)  The  problem  of  opt imal  meshes.  In  the  neighborhood  of  corners. 
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refined  meshes  similar  to  the  ones  discussed  in  the  previous  sections  (radical 
or  geometric)  should  be  used.  In  the  neighborhood  of  edges,  the  use  of 
"needle"  elements  analogous  to  those  used  in  the  neighborhood  of  edges  in 
three  dimensions  is  of  importance  to  resolve  boundary  layers.  Let  us  mention, 
however,  that  the  singularity  here  is  of  different  character  than  the  edge 
singularity  in  three-dimensional  problems.  The  problem  of  the  optimal  mesh 
has  not  been  satisfactory  resolved  as  yet. 

c)  The  pollution  problem.  The  error  in  the  neighborhood  of  the 
boundary  can  now  influence  the  accuracy  of  the  solution  far  from  the  boundary. 
Essentially,  there  is  pollution  present. 

In  addition  to  the  problem  of  plates  discussed  above,  let  us  mention 
that  a  similar  treatment  is  possible  for  other  problems  over  "thin"  domains, 
like  those  arising,  for  instance,  in  shell  theory. 

8.  IMPLEMENTATIONAL  ASPECTS 

In  this  paper,  our  emphasis  has  been  mainly  on  theoretical  aspects  of 
the  p  and  h-p  versions.  In  this  section,  we  briefly  focus  on 
implementational  aspects  which,  of  course,  are  of  essential  importance  in 
terms  of  practical  use. 

As  mentioned  in  the  introduction,  currently  there  are  only  a  few  large 
research  and  commercial  codes  based  on  the  p  and  h-p  versions.  STRIPE 
(Aeronautical  Research  Institute  of  Sweden)  is  a  3-d  adaptive  p  version 
research  code,  Applied  Structure  (Rasna  Corp.  CA,  USA)  and  MSC/PROBE  (MacNeal 
Schwendler,  CA,  USA)  are  p  version  commercial  codes  and  PHLEX  (Computational 
Mechanics,  TX,  USA)  is  an  adaptive  h-p  version  commercial  program.  The 
implementation  of  such  p  and  h-p  version  codes  is  significantly  different  from 
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that  of  h  version  codes.  Let  us  summarize  some  of  the  basic  differences. 

a)  Mesh  construction.  This  is  much  simpler  for  the  p  version,  which 
uses  a  relatively  small  number  of  elements  of  large  size.  With  p  version 
codes,  it  is  important  that  the  user  select  a  good  mesh  to  obtain  a  high  rate 
of  convergence  (See  Remark  4.1). 

b)  Sparsity.  The  local  stiffness  matrices  are  much  less  sparse  for  the 
p  version  than  the  h  version.  Hence,  the  computation  of  stiffness  matrices 
has  to  be  given  special  care,  for  instance  by  taking  into  account  the  computer 
architecture  used.  For  example,  the  program  STRIPE  successfully  exploits  the 
architecture  of  the  CRAY  computer  to  deal  with  this  problem.  Reduced 
sparsity,  especially  in  the  case  of  three-dimensional  problems  with  more  than 
100,000  -  500,000  unknowns,  has  serious  implications  for  the  I/O  time  of 
communication  between  disk  and  fast  memory. 

c)  Adaptivi tv.  In  the  case  of  the  p  version,  the  structure  of  an 
adaptive  program  is  much  simpler  than  the  h  version,  since  only  new  shape 
functions  are  added,  without  changing  the  mesh.  The  architecture  for  the 
adaptive  h-p  version  (e.g.  PHLEX)  is  much  more  complex  (see  [DO]).  The 
adaptive  selection  of  shape  functions  is  usually  determined  by  considerations 
based  on  the  error  in  the  energy.  Essentially,  those  shape  functions  are 
added  which  maximally  change  the  computed  energy.  The  reason  that  adaptive 
principles  for  the  h-p  version  are  more  complex  is  that  simultaneous  decisions 
about  mesh  changes  and  shape  function  selection  have  to  be  made. 

d)  Accuracy  assessment .  The  p  version  allows  the  user  to  control  the 
desired  accuracy  of  the  data  of  interest.  This  control  is  essentially  based 
on  comparison  of  •'he  computed  data  by  increasing  the  degree  p,  so  that  the 
user  can  a  priori  specify  the  required  accuracy.  This  feature  has  been 
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incorporated  e.g.  in  Applied  Structure — see  also  [AB] . 


e)  Linear  solver .  This  can  either  be  direct,  when  many  right  hand 
sides  (loads)  are  present,  or  iterative,  when  the  number  of  unknowns  is  large 
(>  100k,  say)  (see  e.g.  [Mn] ) .  An  effective  preconditioner  can  be  based  on 
the  low  p  (p  =  1,2)  discretization. 

f)  Hierarchic  bases.  As  mentioned  in  Section  2.2,  the  shape  functions 
for  the  p  version  have  a  hierarchic  character.  This  is  exploited  in  adaptive 
computations  and  in  the  construction  of  preconditioners. 

g)  Overal 1  effectivity.  The  effectivity  of  any  methods,  h,  p  or  h-p, 
depends  upon  many  factors.  In  this  paper,  we  have  focused  primarily  on  one 
factor,  the  asymptotic  rate  of  convergence  (with  respect  to  the  number  of 
degrees  of  freedom).  Some  other  important  factors  are  computational  cost, 
man-power  cost  for  data  preparation,  and  proper  usage.  Further,  estimations 
of  efficiency  must  also  take  into  account  various  features  that  may  be 
present,  such  as  adaptive  features,  a  posteriori  error  estimation,  etc.  (The 
capability  of  a  posteriori  error  estimations  is  essential  if  a  program  is  to 
give  reliable  results. ) 

Let  us  mention  with  respect  to  overall  effectivity  that  the  complex 
question  of  comparison  of  different  methods  and  codes  is  usually  accomplished 
through  various  benchmark  problems.  In  this  regard,  an  interesting 
engineering  comparison  study  is  presented  in  ICM] ,  where  the  p  version  is 
compared  to  a  commercial  h  version  code  for  two  test  problems  of  typical 
industrial  complexity.  The  results  presented  underscore  the  viability  of  the 
p  version  not  only  in  comparisons  based  on  the  number  of  degrees  of  freedom, 
but  also  in  those  based  on  the  total  CPU  time. 
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